home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 19 / Mac Magazin and MacEasy Magazine CD - Issue 19.iso / Wissenschaft & Technik / FFTs for RISC 1.1 / fftlib.c < prev    next >
Text File  |  1996-01-30  |  58KB  |  2,433 lines
  1. /* FFT library    */
  2. /* Library of in-place fast fourier transforms    */
  3. /* Forward and inverse complex transforms    */
  4. /* and real forward transform    */
  5. /* Optimized for RISC processors with many registers */
  6. /* Version 1.1 John Green NUWC New London CT    January 96    */
  7. /* Version 1.1 renamed as fftlib from fftbig    */
  8. /* Version 1.1 removed (float *)((char *) ptr) optimization    */
  9. /* Version 1.0  John Green NUWC New London CT    December 95    */
  10. /* (John Green) green_jt@vsdec.nl.nuwc.navy.mil    */
  11. /* green_jt@vsdec.nl.nuwc.navy.mil    */
  12.  
  13. #include <math.h>
  14. #include "fftlib.h"
  15.  
  16. #define MAXMROOT    9    /* 2^(MAXMROOT-1) = # of elements in BRcnt */
  17.  
  18. /* Bit reversed counter */
  19. static const unsigned char BRcnt[256] = {
  20.   0,    128,     64,    192,    32,     160,     96,    224,    16,     144,     80,    208,
  21.  48,    176,    112,    240,     8,     136,     72,    200,    40,     168,    104,    232,
  22.  24,    152,     88,    216,    56,     184,    120,    248,     4,     132,     68,    196,
  23.  36,    164,    100,    228,    20,     148,     84,    212,    52,     180,    116,    244,
  24.  12,    140,     76,    204,    44,     172,    108,    236,    28,     156,     92,    220,
  25.  60,    188,    124,    252,     2,     130,     66,    194,    34,     162,     98,    226,
  26.  18,    146,     82,    210,    50,     178,    114,    242,    10,     138,     74,    202,
  27.  42,    170,    106,    234,    26,     154,     90,    218,    58,     186,    122,    250,
  28.   6,    134,     70,    198,    38,     166,    102,    230,    22,     150,     86,    214,
  29.  54,    182,    118,    246,    14,     142,     78,    206,    46,     174,    110,    238,
  30.  30,    158,     94,    222,    62,     190,    126,    254,     1,     129,     65,    193,
  31.  33,    161,     97,    225,    17,     145,     81,    209,    49,     177,    113,    241,
  32.   9,    137,     73,    201,    41,     169,    105,    233,    25,     153,     89,    217,
  33.  57,    185,    121,    249,     5,     133,     69,    197,    37,     165,    101,    229,
  34.  21,    149,     85,    213,    53,     181,    117,    245,    13,     141,     77,    205,
  35.  45,    173,    109,    237,    29,     157,     93,    221,    61,     189,    125,    253,
  36.   3,    131,     67,    195,    35,     163,     99,    227,    19,     147,     83,    211,
  37.  51,    179,    115,    243,    11,     139,     75,    203,    43,     171,    107,    235,
  38.  27,    155,     91,    219,    59,     187,    123,    251,     7,     135,     71,    199,
  39.  39,    167,    103,    231,    23,     151,     87,    215,    55,     183,    119,    247,
  40.  15,    143,     79,    207,    47,     175,    111,    239,    31,     159,     95,    223,
  41.  63,    191,    127,    255    };
  42.  
  43. /* Table of powers of 2 */
  44. static const long Ntbl[21] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048,
  45.                  4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576};
  46.  
  47. long FFTInit(long *fftMptr, long fftN, float *Utbl){
  48. /* Compute cosine table and check size for complex ffts    */
  49. /* INPUTS */
  50. /* fftN = size of fft    */
  51. /* OUTPUTS */
  52. /* *fftMptr = log2 of fft size    */
  53. /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive)    */
  54. /* RETURNS */
  55. /* 1 if fftN is invalid, 0 otherwise    */
  56. long     i1, ErrVal;
  57. ErrVal = 0;
  58. *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
  59. if ((*fftMptr >= 3) & (*fftMptr <= 19) & (int)(pow(2,*fftMptr)+.5) == fftN)
  60.     for (i1 = 0; i1 <= fftN/4; i1++)
  61.         Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) /  fftN );
  62. else
  63.     ErrVal = 1;
  64. return ErrVal;
  65. }
  66.  
  67. long rFFTInit(long *fftMptr, long fftN, float *Utbl){
  68. /* Compute cosine table and check size for a real input fft    */
  69. /* INPUTS */
  70. /* fftN = size of fft    */
  71. /* OUTPUTS */
  72. /* *fftMptr = log2 of fft size    */
  73. /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive)    */
  74. /* RETURNS */
  75. /* 1 if fftN is invalid, 0 otherwise    */
  76. long     i1, ErrVal;
  77. ErrVal = 0;
  78. *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
  79. if ((*fftMptr >= 4) & (*fftMptr <= 20) & (int)(pow(2,*fftMptr)+.5) == fftN)
  80.     
  81.     for (i1 = 0; i1 <= fftN/4; i1++)
  82.         Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) /  fftN );
  83. else
  84.     ErrVal = 1;
  85. return ErrVal;
  86. }
  87.  
  88. void ffts(float *ioptr, long M, long Rows, float *Utbl){
  89. /* Compute complex fft on the rows of the input array    */
  90. /* INPUTS */
  91. /* M = log2 of fft size    */
  92. /* *ioptr = input data array    */
  93. /* *Utbl = cosine table    */
  94. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  95. /* OUTPUTS */
  96. /* *ioptr = output data array    */
  97.  
  98. long     Flyinc;
  99. long    FlyOffsetA;
  100. long    FlyOffsetAIm;
  101. long    FlyOffsetB;
  102. long    FlyOffsetBIm;
  103. long     NsameU1;
  104. long     NsameU2;
  105. long     NsameU4;
  106. long     diffUcnt;
  107. long     LoopCnt;
  108.  
  109. float    scale;
  110. float    fly0r;
  111. float    fly0i;
  112. float    fly1r;
  113. float    fly1i;
  114. float    fly2r;
  115. float    fly2i;
  116. float    fly3r;
  117. float    fly3i;
  118. float    fly4r;
  119. float    fly4i;
  120. float    fly5r;
  121. float    fly5i;
  122. float    fly6r;
  123. float    fly6i;
  124. float    fly7r;
  125. float    fly7i;
  126. float    U0r;
  127. float    U0i;
  128. float    U1r;
  129. float    U1i;
  130. float    U2r;
  131. float    U2i;
  132. float    U3r;
  133. float    U3i;
  134. float    t0r;
  135. float    t0i;
  136. float    t1r;
  137. float    t1i;
  138.  
  139. float    *fly0P;
  140. float    *fly1P;
  141. float    *fly2P;
  142. float    *fly3P;
  143.  
  144. float    *U0rP;
  145. float    *U0iP;
  146. float    *U1rP;
  147. float    *U1iP;
  148. float    *U2rP;
  149. float    *U2iP;
  150. float    *IOP;
  151. long    U3offset;
  152.  
  153. long     stage;
  154. long     NdiffU;
  155. long     LoopN;
  156.  
  157. const long BRshift = MAXMROOT - (M>>1);
  158. const long Nrems2 = Ntbl[M-(M>>1)+1];
  159. const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)])*2;
  160.  
  161. scale = 2.0;
  162. for (;Rows>0;Rows--){
  163.  
  164. FlyOffsetA = Ntbl[M] * (2/2);
  165. FlyOffsetAIm = FlyOffsetA + 1;
  166. FlyOffsetB = FlyOffsetA + 2;
  167. FlyOffsetBIm = FlyOffsetB + 1;
  168.  
  169. /* BitrevR2 ** bit reverse and first radix 2 stage ******/
  170.  
  171. for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
  172.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  173.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  174.         fly0P = ioptr+ Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
  175.         IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
  176.         fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
  177.         fly0r = *(fly0P);
  178.         fly0i = *(fly0P+1);
  179.         fly1r = *(fly0P+FlyOffsetA);
  180.         fly1i = *(fly0P+FlyOffsetAIm);
  181.         for (; LoopCnt > LoopN;){
  182.             fly2r = *(fly0P+2);
  183.             fly2i = *(fly0P+(2+1));
  184.             fly3r = *(fly0P+FlyOffsetB);
  185.             fly3i = *(fly0P+FlyOffsetBIm);
  186.             fly4r = *(fly1P);
  187.             fly4i = *(fly1P+1);
  188.             fly5r = *(fly1P+FlyOffsetA);
  189.             fly5i = *(fly1P+FlyOffsetAIm);
  190.             fly6r = *(fly1P+2);
  191.             fly6i = *(fly1P+(2+1));
  192.             fly7r = *(fly1P+FlyOffsetB);
  193.             fly7i = *(fly1P+FlyOffsetBIm);
  194.             
  195.             t0r    = fly0r + fly1r;
  196.             t0i    = fly0i + fly1i;
  197.             fly1r = fly0r - fly1r;
  198.             fly1i = fly0i - fly1i;
  199.             t1r = fly2r + fly3r;
  200.             t1i = fly2i + fly3i;
  201.             fly3r = fly2r - fly3r;
  202.             fly3i = fly2i - fly3i;
  203.             fly0r = fly4r + fly5r;
  204.             fly0i = fly4i + fly5i;
  205.             fly5r = fly4r - fly5r;
  206.             fly5i = fly4i - fly5i;
  207.             fly2r = fly6r + fly7r;
  208.             fly2i = fly6i + fly7i;
  209.             fly7r = fly6r - fly7r;
  210.             fly7i = fly6i - fly7i;
  211.  
  212.             *(fly1P) = t0r;
  213.             *(fly1P+1) = t0i;
  214.             *(fly1P+2) = fly1r;
  215.             *(fly1P+(2+1)) = fly1i;
  216.             *(fly1P+FlyOffsetA) = t1r;
  217.             *(fly1P+FlyOffsetAIm) = t1i;
  218.             *(fly1P+FlyOffsetB) = fly3r;
  219.             *(fly1P+FlyOffsetBIm) = fly3i;
  220.             *(fly0P) = fly0r;
  221.             *(fly0P+1) = fly0i;
  222.             *(fly0P+2) = fly5r;
  223.             *(fly0P+(2+1)) = fly5i;
  224.             *(fly0P+FlyOffsetA) = fly2r;
  225.             *(fly0P+FlyOffsetAIm) = fly2i;
  226.             *(fly0P+FlyOffsetB) = fly7r;
  227.             *(fly0P+FlyOffsetBIm) = fly7i;
  228.  
  229.             fly0P -= Nrems2;
  230.             fly0r = *(fly0P);
  231.             fly0i = *(fly0P+1);
  232.             fly1r = *(fly0P+FlyOffsetA);
  233.             fly1i = *(fly0P+FlyOffsetAIm);
  234.             LoopCnt -= 1;
  235.             fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
  236.         };
  237.         fly2r = *(fly0P+2);
  238.         fly2i = *(fly0P+(2+1));
  239.         fly3r = *(fly0P+FlyOffsetB);
  240.         fly3i = *(fly0P+FlyOffsetBIm);
  241.  
  242.         t0r   = fly0r + fly1r;
  243.         t0i   = fly0i + fly1i;
  244.         fly1r = fly0r - fly1r;
  245.         fly1i = fly0i - fly1i;
  246.         t1r   = fly2r + fly3r;
  247.         t1i   = fly2i + fly3i;
  248.         fly3r = fly2r - fly3r;
  249.         fly3i = fly2i - fly3i;
  250.  
  251.         *(fly0P) = t0r;
  252.         *(fly0P+1) = t0i;
  253.         *(fly0P+2) = fly1r;
  254.         *(fly0P+(2+1)) = fly1i;
  255.         *(fly0P+FlyOffsetA) = t1r;
  256.         *(fly0P+FlyOffsetAIm) = t1i;
  257.         *(fly0P+FlyOffsetB) = fly3r;
  258.         *(fly0P+FlyOffsetBIm) = fly3i;
  259.  
  260.     };
  261. };
  262.  
  263.  
  264.  
  265. /**** FFTC  **************/
  266.  
  267. NdiffU = 2;
  268. Flyinc =  (NdiffU);
  269.  
  270. NsameU4 = Ntbl[M]/4;
  271. LoopN = Ntbl[M-3];
  272.  
  273. stage = ((M-1)/3);
  274. if ( (M-1-(stage * 3)) != 0 ){
  275.     FlyOffsetA =  Flyinc << 2;
  276.     FlyOffsetB =  FlyOffsetA << 1;
  277.     FlyOffsetAIm =  FlyOffsetA + 1;
  278.     FlyOffsetBIm =  FlyOffsetB + 1;
  279.     if ( (M-1-(stage * 3)) == 1 ){
  280.         /* 1 radix 2 stage */
  281.  
  282.         IOP = ioptr;
  283.         fly0P = IOP;
  284.         fly1P = (IOP+Flyinc);
  285.         fly2P = (fly1P+Flyinc);
  286.         fly3P = (fly2P+Flyinc);
  287.         
  288.             /* Butterflys        */
  289.             /*
  290.             t0    -    -    t0
  291.             t1    -    -    t1
  292.             f2    -  1-    f5
  293.             f1    - -i-    f7
  294.             f4    -    -    f4
  295.             f0    -    -    f0
  296.             f6    -  1-    f2
  297.             f3    - -i-    f1
  298.             */
  299.  
  300.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  301.             t0r = *(fly0P);
  302.             t0i = *(fly0P + 1);
  303.             t1r = *(fly1P);
  304.             t1i = *(fly1P + 1);
  305.             fly2r = *(fly2P);
  306.             fly2i = *(fly2P + 1);
  307.             fly1r = *(fly3P);
  308.             fly1i = *(fly3P + 1);
  309.             fly4r = *(fly0P + FlyOffsetA);
  310.             fly4i = *(fly0P + FlyOffsetAIm);
  311.             fly0r = *(fly1P + FlyOffsetA);
  312.             fly0i = *(fly1P + FlyOffsetAIm);
  313.             fly6r = *(fly2P + FlyOffsetA);
  314.             fly6i = *(fly2P + FlyOffsetAIm);
  315.             fly3r = *(fly3P + FlyOffsetA);
  316.             fly3i = *(fly3P + FlyOffsetAIm);
  317.         
  318.             fly5r = t0r - fly2r;
  319.             fly5i = t0i - fly2i;
  320.             t0r = t0r + fly2r;
  321.             t0i = t0i + fly2i;
  322.  
  323.             fly7r = t1r - fly1i;
  324.             fly7i = t1i + fly1r;
  325.             t1r = t1r + fly1i;
  326.             t1i = t1i - fly1r;
  327.  
  328.             fly2r = fly4r - fly6r;
  329.             fly2i = fly4i - fly6i;
  330.             fly4r = fly4r + fly6r;
  331.             fly4i = fly4i + fly6i;
  332.  
  333.             fly1r = fly0r - fly3i;
  334.             fly1i = fly0i + fly3r;
  335.             fly0r = fly0r + fly3i;
  336.             fly0i = fly0i - fly3r;
  337.  
  338.             *(fly2P) = fly5r;
  339.             *(fly2P + 1) = fly5i;
  340.             *(fly0P) = t0r;
  341.             *(fly0P + 1) = t0i;
  342.             *(fly3P) = fly7r;
  343.             *(fly3P + 1) = fly7i;
  344.             *(fly1P) = t1r;
  345.             *(fly1P + 1) = t1i;
  346.             *(fly2P + FlyOffsetA) = fly2r;
  347.             *(fly2P + FlyOffsetAIm) = fly2i;
  348.             *(fly0P + FlyOffsetA) = fly4r;
  349.             *(fly0P + FlyOffsetAIm) = fly4i;
  350.             *(fly3P + FlyOffsetA) = fly1r;
  351.             *(fly3P + FlyOffsetAIm) = fly1i;
  352.             *(fly1P + FlyOffsetA) = fly0r;
  353.             *(fly1P + FlyOffsetAIm) = fly0i;
  354.  
  355.             fly0P = (fly0P + FlyOffsetB);
  356.             fly1P = (fly1P + FlyOffsetB);
  357.             fly2P = (fly2P + FlyOffsetB);
  358.             fly3P = (fly3P + FlyOffsetB);
  359.         };
  360.  
  361.         NsameU4 >>= 1;
  362.         LoopN >>= 1;
  363.         NdiffU <<= 1;
  364.         Flyinc = Flyinc << 1;
  365.     }
  366.     else{
  367.         /* 1 radix 4 stage */
  368.         IOP = ioptr;
  369.  
  370.         U3r =  0.7071067811865475244008443621; /* sqrt(0.5);    */    
  371.         U3i = U3r;    
  372.         fly0P = IOP;
  373.         fly1P = (IOP+Flyinc);
  374.         fly2P = (fly1P+Flyinc);
  375.         fly3P = (fly2P+Flyinc);
  376.         
  377.             /* Butterflys        */
  378.             /*
  379.             t0    -    -    t0    -    -    t0
  380.             t1    -    -    t1    -    -    t1
  381.             f2    -  1-    f5    -    -    f5
  382.             f1    - -i-    f7    -    -    f7
  383.             f4    -    -    f4    -  1-    f6
  384.             f0    -    -    f0    -U3    -    f3
  385.             f6    -  1-    f2    - -i-    f4
  386.             f3    - -i-    f1    -U3a-    f2
  387.             */
  388.  
  389.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  390.             t0r = *(fly0P);
  391.             t0i = *(fly0P + 1);
  392.             t1r = *(fly1P);
  393.             t1i = *(fly1P + 1);
  394.             fly2r = *(fly2P);
  395.             fly2i = *(fly2P + 1);
  396.             fly1r = *(fly3P);
  397.             fly1i = *(fly3P + 1);
  398.             fly4r = *(fly0P + FlyOffsetA);
  399.             fly4i = *(fly0P + FlyOffsetAIm);
  400.             fly0r = *(fly1P + FlyOffsetA);
  401.             fly0i = *(fly1P + FlyOffsetAIm);
  402.             fly6r = *(fly2P + FlyOffsetA);
  403.             fly6i = *(fly2P + FlyOffsetAIm);
  404.             fly3r = *(fly3P + FlyOffsetA);
  405.             fly3i = *(fly3P + FlyOffsetAIm);
  406.     
  407.             fly5r = t0r - fly2r;
  408.             fly5i = t0i - fly2i;
  409.             t0r = t0r + fly2r;
  410.             t0i = t0i + fly2i;
  411.     
  412.             fly7r = t1r - fly1i;
  413.             fly7i = t1i + fly1r;
  414.             t1r = t1r + fly1i;
  415.             t1i = t1i - fly1r;
  416.     
  417.             fly2r = fly4r - fly6r;
  418.             fly2i = fly4i - fly6i;
  419.             fly4r = fly4r + fly6r;
  420.             fly4i = fly4i + fly6i;
  421.     
  422.             fly1r = fly0r - fly3i;
  423.             fly1i = fly0i + fly3r;
  424.             fly0r = fly0r + fly3i;
  425.             fly0i = fly0i - fly3r;
  426.     
  427.     
  428.             fly6r = t0r - fly4r;
  429.             fly6i = t0i - fly4i;
  430.             t0r = t0r + fly4r;
  431.             t0i = t0i + fly4i;
  432.     
  433.             fly3r = fly5r - fly2i;
  434.             fly3i = fly5i + fly2r;
  435.             fly5r = fly5r + fly2i;
  436.             fly5i = fly5i - fly2r;
  437.     
  438.             fly4r = t1r - U3r * fly0r;
  439.             fly4r = fly4r - U3i * fly0i;
  440.             fly4i = t1i + U3i * fly0r;
  441.             fly4i = fly4i - U3r * fly0i;
  442.             t1r = scale * t1r - fly4r;
  443.             t1i = scale * t1i - fly4i;
  444.     
  445.             fly2r = fly7r + U3i * fly1r;
  446.             fly2r = fly2r - U3r * fly1i;
  447.             fly2i = fly7i + U3r * fly1r;
  448.             fly2i = fly2i + U3i * fly1i;
  449.             fly7r = scale * fly7r - fly2r;
  450.             fly7i = scale * fly7i - fly2i;
  451.     
  452.             *(fly0P + FlyOffsetA) = fly6r;
  453.             *(fly0P + FlyOffsetAIm) = fly6i;
  454.             *(fly0P) = t0r;
  455.             *(fly0P + 1) = t0i;
  456.             *(fly2P + FlyOffsetA) = fly3r;
  457.             *(fly2P + FlyOffsetAIm) = fly3i;
  458.             *(fly2P) = fly5r;
  459.             *(fly2P + 1) = fly5i;
  460.             *(fly1P + FlyOffsetA) = fly4r;
  461.             *(fly1P + FlyOffsetAIm) = fly4i;
  462.             *(fly1P) = t1r;
  463.             *(fly1P + 1) = t1i;
  464.             *(fly3P + FlyOffsetA) = fly2r;
  465.             *(fly3P + FlyOffsetAIm) = fly2i;
  466.             *(fly3P) = fly7r;
  467.             *(fly3P + 1) = fly7i;
  468.         
  469.             fly0P = (fly0P + FlyOffsetB);
  470.             fly1P = (fly1P + FlyOffsetB);
  471.             fly2P = (fly2P + FlyOffsetB);
  472.             fly3P = (fly3P + FlyOffsetB);
  473.  
  474.         };
  475.  
  476.         NsameU4 >>= 2;
  477.         LoopN >>= 2;
  478.         NdiffU <<= 2;
  479.         Flyinc = Flyinc << 2;
  480.     };
  481. };
  482.  
  483. NsameU2 = NsameU4 >> 1;
  484. NsameU1 = NsameU2 >> 1;
  485. Flyinc <<= 1;
  486. FlyOffsetA =  Flyinc << 2;
  487. FlyOffsetB =  FlyOffsetA << 1;
  488. FlyOffsetAIm =  FlyOffsetA + 1;
  489. FlyOffsetBIm =  FlyOffsetB + 1;
  490. LoopN >>= 1;
  491. /*   ****** RADIX 8 Stages    */
  492. for (stage = stage<<1; stage > 0 ; stage--){
  493.  
  494.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  495.  
  496. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  497.     if(!(stage&1)){
  498.         U0rP = &Utbl[0];
  499.         U0iP = &Utbl[Ntbl[M-2]];
  500.         U1rP = U0rP;
  501.         U1iP = U0iP;
  502.         U2rP = U0rP;
  503.         U2iP = U0iP;
  504.         U3offset = (Ntbl[M]) / 8;
  505.  
  506.         IOP = ioptr;
  507.     
  508.         U0r =  *U0rP,
  509.         U0i =  *U0iP;
  510.         U1r =  *U1rP,
  511.         U1i =  *U1iP;
  512.         U2r =  *U2rP,
  513.         U2i =  *U2iP;
  514.         U3r =  *( U2rP + U3offset);
  515.         U3i =  *( U2iP - U3offset);
  516.     }
  517.     
  518.     fly0P = IOP;
  519.     fly1P = (IOP+Flyinc);
  520.     fly2P = (fly1P+Flyinc);
  521.     fly3P = (fly2P+Flyinc);
  522.  
  523.     for (diffUcnt = (NdiffU)>>1; diffUcnt != 0; diffUcnt--){
  524.  
  525.             /* Butterflys        */
  526.             /*
  527.             f0    -    -    t0    -    -    t0    -    -    t0
  528.             f1    -U0    -    t1    -    -    t1    -    -    t1
  529.             f2    -    -    f2    -U1    -    f5    -    -    f3
  530.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  531.             f4    -    -    f4    -    -    f4    -U2    -    f6
  532.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  533.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  534.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  535.             */
  536.         
  537.         fly0r = *(IOP);
  538.         fly0i = *(IOP+1);
  539.         fly1r = *(fly1P);
  540.         fly1i = *(fly1P+1);
  541.  
  542.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  543.     
  544.             fly2r = *(fly2P);
  545.             fly2i = *(fly2P + 1);
  546.             fly3r = *(fly3P);
  547.             fly3i = *(fly3P + 1);
  548.             fly4r = *(fly0P + FlyOffsetA);
  549.             fly4i = *(fly0P + FlyOffsetAIm);
  550.             fly5r = *(fly1P + FlyOffsetA);
  551.             fly5i = *(fly1P + FlyOffsetAIm);
  552.             fly6r = *(fly2P + FlyOffsetA);
  553.             fly6i = *(fly2P + FlyOffsetAIm);
  554.             fly7r = *(fly3P + FlyOffsetA);
  555.             fly7i = *(fly3P + FlyOffsetAIm);
  556.  
  557.             t1r = fly0r - U0r * fly1r;
  558.             t1r = t1r - U0i * fly1i;
  559.             t1i = fly0i + U0i * fly1r;
  560.             t1i = t1i - U0r * fly1i;
  561.             t0r = scale * fly0r - t1r;
  562.             t0i = scale * fly0i - t1i;
  563.     
  564.             fly1r = fly2r - U0r * fly3r;
  565.             fly1r = fly1r - U0i * fly3i;
  566.             fly1i = fly2i + U0i * fly3r;
  567.             fly1i = fly1i - U0r * fly3i;
  568.             fly2r = scale * fly2r - fly1r;
  569.             fly2i = scale * fly2i - fly1i;
  570.     
  571.             fly0r = fly4r - U0r * fly5r;
  572.             fly0r = fly0r - U0i * fly5i;
  573.             fly0i = fly4i + U0i * fly5r;
  574.             fly0i = fly0i - U0r * fly5i;
  575.             fly4r = scale * fly4r - fly0r;
  576.             fly4i = scale * fly4i - fly0i;
  577.     
  578.             fly3r = fly6r - U0r * fly7r;
  579.             fly3r = fly3r - U0i * fly7i;
  580.             fly3i = fly6i + U0i * fly7r;
  581.             fly3i = fly3i - U0r * fly7i;
  582.             fly6r = scale * fly6r - fly3r;
  583.             fly6i = scale * fly6i - fly3i;
  584.     
  585.  
  586.             fly5r = t0r - U1r * fly2r;
  587.             fly5r = fly5r - U1i * fly2i;
  588.             fly5i = t0i + U1i * fly2r;
  589.             fly5i = fly5i - U1r * fly2i;
  590.             t0r = scale * t0r - fly5r;
  591.             t0i = scale * t0i - fly5i;
  592.  
  593.             fly7r = t1r + U1i * fly1r;
  594.             fly7r = fly7r - U1r * fly1i;
  595.             fly7i = t1i + U1r * fly1r;
  596.             fly7i = fly7i + U1i * fly1i;
  597.             t1r = scale * t1r - fly7r;
  598.             t1i = scale * t1i - fly7i;
  599.  
  600.             fly2r = fly4r - U1r * fly6r;
  601.             fly2r = fly2r - U1i * fly6i;
  602.             fly2i = fly4i + U1i * fly6r;
  603.             fly2i = fly2i - U1r * fly6i;
  604.             fly4r = scale * fly4r - fly2r;
  605.             fly4i = scale * fly4i - fly2i;
  606.  
  607.             fly1r = fly0r + U1i * fly3r;
  608.             fly1r = fly1r - U1r * fly3i;
  609.             fly1i = fly0i + U1r * fly3r;
  610.             fly1i = fly1i + U1i * fly3i;
  611.             fly0r = scale * fly0r - fly1r;
  612.             fly0i = scale * fly0i - fly1i;
  613.  
  614.             fly6r = t0r - U2r * fly4r;
  615.             fly6r = fly6r - U2i * fly4i;
  616.             fly6i = t0i + U2i * fly4r;
  617.             fly6i = fly6i - U2r * fly4i;
  618.             t0r = scale * t0r - fly6r;
  619.             t0i = scale * t0i - fly6i;
  620.  
  621.             fly3r = fly5r - U2i * fly2r;
  622.             fly3r = fly3r + U2r * fly2i;
  623.             fly3i = fly5i - U2r * fly2r;
  624.             fly3i = fly3i - U2i * fly2i;
  625.             fly2r = scale * fly5r - fly3r;
  626.             fly2i = scale * fly5i - fly3i;
  627.  
  628.             fly4r = t1r - U3r * fly0r;
  629.             fly4r = fly4r - U3i * fly0i;
  630.             fly4i = t1i + U3i * fly0r;
  631.             fly4i = fly4i - U3r * fly0i;
  632.             t1r = scale * t1r - fly4r;
  633.             t1i = scale * t1i - fly4i;
  634.  
  635.             fly5r = fly7r + U3i * fly1r;
  636.             fly5r = fly5r - U3r * fly1i;
  637.             fly5i = fly7i + U3r * fly1r;
  638.             fly5i = fly5i + U3i * fly1i;
  639.             fly7r = scale * fly7r - fly5r;
  640.             fly7i = scale * fly7i - fly5i;
  641.  
  642.             *(fly0P + FlyOffsetA) = fly6r;
  643.             *(fly0P + FlyOffsetAIm) = fly6i;
  644.             *(fly0P) = t0r;
  645.             *(fly0P + 1) = t0i;
  646.             *(fly2P) = fly3r;
  647.             *(fly2P + 1) = fly3i;
  648.             *(fly2P + FlyOffsetA) = fly2r;
  649.             *(fly2P + FlyOffsetAIm) = fly2i;
  650.  
  651.             fly0r = *(fly0P + FlyOffsetB);
  652.             fly0i = *(fly0P + FlyOffsetBIm);
  653.  
  654.             *(fly1P + FlyOffsetA) = fly4r;
  655.             *(fly1P + FlyOffsetAIm) = fly4i;
  656.             *(fly1P) = t1r;
  657.             *(fly1P + 1) = t1i;
  658.  
  659.             fly1r = *(fly1P + FlyOffsetB);
  660.             fly1i = *(fly1P + FlyOffsetBIm);
  661.  
  662.             *(fly3P + FlyOffsetA) = fly5r;
  663.             *(fly3P + FlyOffsetAIm) = fly5i;
  664.             *(fly3P) = fly7r;
  665.             *(fly3P + 1) = fly7i;
  666.  
  667.             fly0P = (fly0P + FlyOffsetB);
  668.             fly1P = (fly1P + FlyOffsetB);
  669.             fly2P = (fly2P + FlyOffsetB);
  670.             fly3P = (fly3P + FlyOffsetB);
  671.         };
  672.         fly2r = *(fly2P);
  673.         fly2i = *(fly2P + 1);
  674.         fly3r = *(fly3P);
  675.         fly3i = *(fly3P + 1);
  676.         fly4r = *(fly0P + FlyOffsetA);
  677.         fly4i = *(fly0P + FlyOffsetAIm);
  678.         fly5r = *(fly1P + FlyOffsetA);
  679.         fly5i = *(fly1P + FlyOffsetAIm);
  680.         fly6r = *(fly2P + FlyOffsetA);
  681.         fly6i = *(fly2P + FlyOffsetAIm);
  682.         fly7r = *(fly3P + FlyOffsetA);
  683.         fly7i = *(fly3P + FlyOffsetAIm);
  684.  
  685.         t1r = fly0r - U0r * fly1r;
  686.         t1r = t1r - U0i * fly1i;
  687.         t1i = fly0i + U0i * fly1r;
  688.         t1i = t1i - U0r * fly1i;
  689.         t0r = scale * fly0r - t1r;
  690.         t0i = scale * fly0i - t1i;
  691.  
  692.         fly1r = fly2r - U0r * fly3r;
  693.         fly1r = fly1r - U0i * fly3i;
  694.         fly1i = fly2i + U0i * fly3r;
  695.         fly1i = fly1i - U0r * fly3i;
  696.         fly2r = scale * fly2r - fly1r;
  697.         fly2i = scale * fly2i - fly1i;
  698.  
  699.         fly0r = fly4r - U0r * fly5r;
  700.         fly0r = fly0r - U0i * fly5i;
  701.         fly0i = fly4i + U0i * fly5r;
  702.         fly0i = fly0i - U0r * fly5i;
  703.         fly4r = scale * fly4r - fly0r;
  704.         fly4i = scale * fly4i - fly0i;
  705.  
  706.         fly3r = fly6r - U0r * fly7r;
  707.         fly3r = fly3r - U0i * fly7i;
  708.         fly3i = fly6i + U0i * fly7r;
  709.         fly3i = fly3i - U0r * fly7i;
  710.         fly6r = scale * fly6r - fly3r;
  711.         fly6i = scale * fly6i - fly3i;
  712.  
  713.         fly5r = t0r - U1r * fly2r;
  714.         fly5r = fly5r - U1i * fly2i;
  715.         fly5i = t0i + U1i * fly2r;
  716.         fly5i = fly5i - U1r * fly2i;
  717.         t0r = scale * t0r - fly5r;
  718.         t0i = scale * t0i - fly5i;
  719.  
  720.         fly7r = t1r + U1i * fly1r;
  721.         fly7r = fly7r - U1r * fly1i;
  722.         fly7i = t1i + U1r * fly1r;
  723.         fly7i = fly7i + U1i * fly1i;
  724.         t1r = scale * t1r - fly7r;
  725.         t1i = scale * t1i - fly7i;
  726.  
  727.         fly2r = fly4r - U1r * fly6r;
  728.         fly2r = fly2r - U1i * fly6i;
  729.         fly2i = fly4i + U1i * fly6r;
  730.         fly2i = fly2i - U1r * fly6i;
  731.         fly4r = scale * fly4r - fly2r;
  732.         fly4i = scale * fly4i - fly2i;
  733.  
  734.         fly1r = fly0r + U1i * fly3r;
  735.         fly1r = fly1r - U1r * fly3i;
  736.         fly1i = fly0i + U1r * fly3r;
  737.         fly1i = fly1i + U1i * fly3i;
  738.         fly0r = scale * fly0r - fly1r;
  739.         fly0i = scale * fly0i - fly1i;
  740.  
  741.         U0i = *(U0iP = (U0iP - NsameU4));
  742.         U0r = *(U0rP = (U0rP + NsameU4));        
  743.         U1r = *(U1rP = (U1rP + NsameU2));
  744.         U1i = *(U1iP = (U1iP - NsameU2));
  745.         if(stage&1)
  746.             U0r = -U0r;
  747.  
  748.         fly6r = t0r - U2r * fly4r;
  749.         fly6r = fly6r - U2i * fly4i;
  750.         fly6i = t0i + U2i * fly4r;
  751.         fly6i = fly6i - U2r * fly4i;
  752.         t0r = scale * t0r - fly6r;
  753.         t0i = scale * t0i - fly6i;
  754.  
  755.         fly3r = fly5r - U2i * fly2r;
  756.         fly3r = fly3r + U2r * fly2i;
  757.         fly3i = fly5i - U2r * fly2r;
  758.         fly3i = fly3i - U2i * fly2i;
  759.         fly2r = scale * fly5r - fly3r;
  760.         fly2i = scale * fly5i - fly3i;
  761.  
  762.         fly4r = t1r - U3r * fly0r;
  763.         fly4r = fly4r - U3i * fly0i;
  764.         fly4i = t1i + U3i * fly0r;
  765.         fly4i = fly4i - U3r * fly0i;
  766.         t1r = scale * t1r - fly4r;
  767.         t1i = scale * t1i - fly4i;
  768.  
  769.         fly5r = fly7r + U3i * fly1r;
  770.         fly5r = fly5r - U3r * fly1i;
  771.         fly5i = fly7i + U3r * fly1r;
  772.         fly5i = fly5i + U3i * fly1i;
  773.         fly7r = scale * fly7r - fly5r;
  774.         fly7i = scale * fly7i - fly5i;
  775.  
  776.         *(fly0P + FlyOffsetA) = fly6r;
  777.         *(fly0P + FlyOffsetAIm) = fly6i;
  778.         *(fly0P) = t0r;
  779.         *(fly0P + 1) = t0i;
  780.  
  781.         U2r = *(U2rP = (U2rP + NsameU1));
  782.         U2i = *(U2iP = (U2iP - NsameU1));
  783.  
  784.         *(fly2P) = fly3r;
  785.         *(fly2P + 1) = fly3i;
  786.         *(fly2P + FlyOffsetA) = fly2r;
  787.         *(fly2P + FlyOffsetAIm) = fly2i;
  788.         *(fly1P + FlyOffsetA) = fly4r;
  789.         *(fly1P + FlyOffsetAIm) = fly4i;
  790.         *(fly1P) = t1r;
  791.         *(fly1P + 1) = t1i;
  792.  
  793.         U3r = *( U2rP + U3offset);
  794.         U3i = *( U2iP - U3offset);
  795.  
  796.         *(fly3P + FlyOffsetA) = fly5r;
  797.         *(fly3P + FlyOffsetAIm) = fly5i;
  798.         *(fly3P) = fly7r;
  799.         *(fly3P + 1) = fly7i;
  800.  
  801.         IOP = IOP + 2;
  802.         fly0P = IOP;
  803.         fly1P = (IOP+Flyinc);
  804.         fly2P = (fly1P+Flyinc);
  805.         fly3P = (fly2P+Flyinc);
  806.     };
  807.     NsameU4 = -NsameU4;
  808.  
  809.     if(stage&1){
  810.         LoopN >>= 3;
  811.         NsameU1 >>= 3;
  812.         NsameU2 >>= 3;
  813.         NsameU4 >>= 3;
  814.         NdiffU <<= 3;
  815.         Flyinc <<= 3;
  816.         FlyOffsetA <<= 3;
  817.         FlyOffsetB <<= 3;
  818.         FlyOffsetAIm =  FlyOffsetA + 1;
  819.         FlyOffsetBIm =  FlyOffsetB + 1;
  820.     }
  821. }
  822. ioptr += 2*Ntbl[M];
  823. }
  824. }
  825.  
  826. void iffts(float *ioptr, long M, long Rows, float *Utbl){
  827. /* Compute inverse complex fft on the rows of the input array    */
  828. /* INPUTS */
  829. /* M = log2 of fft size    */
  830. /* *ioptr = input data array    */
  831. /* *Utbl = cosine table    */
  832. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  833. /* OUTPUTS */
  834. /* *ioptr = output data array    */
  835.  
  836. long     Flyinc;
  837. long    FlyOffsetA;
  838. long    FlyOffsetAIm;
  839. long    FlyOffsetB;
  840. long    FlyOffsetBIm;
  841. long     NsameU1;
  842. long     NsameU2;
  843. long     NsameU4;
  844. long     diffUcnt;
  845. long     LoopCnt;
  846.  
  847. float    scale;
  848. float    fly0r;
  849. float    fly0i;
  850. float    fly1r;
  851. float    fly1i;
  852. float    fly2r;
  853. float    fly2i;
  854. float    fly3r;
  855. float    fly3i;
  856. float    fly4r;
  857. float    fly4i;
  858. float    fly5r;
  859. float    fly5i;
  860. float    fly6r;
  861. float    fly6i;
  862. float    fly7r;
  863. float    fly7i;
  864. float    U0r;
  865. float    U0i;
  866. float    U1r;
  867. float    U1i;
  868. float    U2r;
  869. float    U2i;
  870. float    U3r;
  871. float    U3i;
  872. float    t0r;
  873. float    t0i;
  874. float    t1r;
  875. float    t1i;
  876.  
  877. float    *fly0P;
  878. float    *fly1P;
  879. float    *fly2P;
  880. float    *fly3P;
  881. float    *U0rP;
  882. float    *U0iP;
  883. float    *U1rP;
  884. float    *U1iP;
  885. float    *U2rP;
  886. float    *U2iP;
  887. float    *IOP;
  888. long    U3offset;
  889.  
  890. long     stage;
  891. long     NdiffU;
  892. long     LoopN;
  893.  
  894. const long BRshift = MAXMROOT - (M>>1);
  895. const long Nrems2 = Ntbl[M-(M>>1)+1];
  896. const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)])*2;
  897.  
  898. for (;Rows>0;Rows--){
  899.  
  900. FlyOffsetA = Ntbl[M] * 2/2;
  901. FlyOffsetAIm = FlyOffsetA + 1;
  902. FlyOffsetB = FlyOffsetA + 2;
  903. FlyOffsetBIm = FlyOffsetB + 1;
  904.  
  905. /* BitrevR2 ** bit reverse and first radix 2 stage ******/
  906.  
  907. scale = 1./Ntbl[M];
  908. for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
  909.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  910.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  911.         fly0P = ioptr + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
  912.         IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
  913.         fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
  914.         fly0r = *(fly0P);
  915.         fly0i = *(fly0P+1);
  916.         fly1r = *(fly0P+FlyOffsetA);
  917.         fly1i = *(fly0P+FlyOffsetAIm);
  918.         for (; LoopCnt > LoopN;){
  919.             fly2r = *(fly0P+2);
  920.             fly2i = *(fly0P+(2+1));
  921.             fly3r = *(fly0P+FlyOffsetB);
  922.             fly3i = *(fly0P+FlyOffsetBIm);
  923.             fly4r = *(fly1P);
  924.             fly4i = *(fly1P+1);
  925.             fly5r = *(fly1P+FlyOffsetA);
  926.             fly5i = *(fly1P+FlyOffsetAIm);
  927.             fly6r = *(fly1P+2);
  928.             fly6i = *(fly1P+(2+1));
  929.             fly7r = *(fly1P+FlyOffsetB);
  930.             fly7i = *(fly1P+FlyOffsetBIm);
  931.             
  932.             t0r   = fly0r + fly1r;
  933.             t0i   = fly0i + fly1i;
  934.             fly1r = fly0r - fly1r;
  935.             fly1i = fly0i - fly1i;
  936.             t1r   = fly2r + fly3r;
  937.             t1i   = fly2i + fly3i;
  938.             fly3r = fly2r - fly3r;
  939.             fly3i = fly2i - fly3i;
  940.             fly0r = fly4r + fly5r;
  941.             fly0i = fly4i + fly5i;
  942.             fly5r = fly4r - fly5r;
  943.             fly5i = fly4i - fly5i;
  944.             fly2r = fly6r + fly7r;
  945.             fly2i = fly6i + fly7i;
  946.             fly7r = fly6r - fly7r;
  947.             fly7i = fly6i - fly7i;
  948.  
  949.             *(fly1P) = scale*t0r;
  950.             *(fly1P+1) = scale*t0i;
  951.             *(fly1P+2) = scale*fly1r;
  952.             *(fly1P+(2+1)) = scale*fly1i;
  953.             *(fly1P+FlyOffsetA) = scale*t1r;
  954.             *(fly1P+FlyOffsetAIm) = scale*t1i;
  955.             *(fly1P+FlyOffsetB) = scale*fly3r;
  956.             *(fly1P+FlyOffsetBIm) = scale*fly3i;
  957.             *(fly0P) = scale*fly0r;
  958.             *(fly0P+1) = scale*fly0i;
  959.             *(fly0P+2) = scale*fly5r;
  960.             *(fly0P+(2+1)) = scale*fly5i;
  961.             *(fly0P+FlyOffsetA) = scale*fly2r;
  962.             *(fly0P+FlyOffsetAIm) = scale*fly2i;
  963.             *(fly0P+FlyOffsetB) = scale*fly7r;
  964.             *(fly0P+FlyOffsetBIm) = scale*fly7i;
  965.  
  966.             fly0P -= Nrems2;
  967.             fly0r = *(fly0P);
  968.             fly0i = *(fly0P+1);
  969.             fly1r = *(fly0P+FlyOffsetA);
  970.             fly1i = *(fly0P+FlyOffsetAIm);
  971.             LoopCnt -= 1;
  972.             fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
  973.         };
  974.         fly2r = *(fly0P+2);
  975.         fly2i = *(fly0P+(2+1));
  976.         fly3r = *(fly0P+FlyOffsetB);
  977.         fly3i = *(fly0P+FlyOffsetBIm);
  978.  
  979.         t0r   = fly0r + fly1r;
  980.         t0i   = fly0i + fly1i;
  981.         fly1r = fly0r - fly1r;
  982.         fly1i = fly0i - fly1i;
  983.         t1r   = fly2r + fly3r;
  984.         t1i   = fly2i + fly3i;
  985.         fly3r = fly2r - fly3r;
  986.         fly3i = fly2i - fly3i;
  987.  
  988.         *(fly0P) = scale*t0r;
  989.         *(fly0P+1) = scale*t0i;
  990.         *(fly0P+2) = scale*fly1r;
  991.         *(fly0P+(2+1)) = scale*fly1i;
  992.         *(fly0P+FlyOffsetA) = scale*t1r;
  993.         *(fly0P+FlyOffsetAIm) = scale*t1i;
  994.         *(fly0P+FlyOffsetB) = scale*fly3r;
  995.         *(fly0P+FlyOffsetBIm) = scale*fly3i;
  996.  
  997.     };
  998. };
  999.  
  1000. /**** FFTC  **************/
  1001.  
  1002. scale = 2.0;
  1003.  
  1004. NdiffU = 2;
  1005. Flyinc =  (NdiffU);
  1006.  
  1007. NsameU4 = Ntbl[M]/4;
  1008. LoopN = Ntbl[M-3];
  1009.  
  1010. stage = ((M-1)/3);
  1011. if ( (M-1-(stage * 3)) != 0 ){
  1012.     FlyOffsetA =  Flyinc << 2;
  1013.     FlyOffsetB =  FlyOffsetA << 1;
  1014.     FlyOffsetAIm =  FlyOffsetA + 1;
  1015.     FlyOffsetBIm =  FlyOffsetB + 1;
  1016.     if ( (M-1-(stage * 3)) == 1 ){
  1017.         /* 1 radix 2 stage */
  1018.  
  1019.         IOP = ioptr;
  1020.         fly0P = IOP;
  1021.         fly1P = (IOP+Flyinc);
  1022.         fly2P = (fly1P+Flyinc);
  1023.         fly3P = (fly2P+Flyinc);
  1024.         
  1025.             /* Butterflys        */
  1026.             /*
  1027.             t0    -    -    t0
  1028.             t1    -    -    t1
  1029.             f2    -  1-    f5
  1030.             f1    - -i-    f7
  1031.             f4    -    -    f4
  1032.             f0    -    -    f0
  1033.             f6    -  1-    f2
  1034.             f3    - -i-    f1
  1035.             */
  1036.  
  1037.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1038.             t0r = *(fly0P);
  1039.             t0i = *(fly0P + 1);
  1040.             t1r = *(fly1P);
  1041.             t1i = *(fly1P + 1);
  1042.             fly2r = *(fly2P);
  1043.             fly2i = *(fly2P + 1);
  1044.             fly1r = *(fly3P);
  1045.             fly1i = *(fly3P + 1);
  1046.             fly4r = *(fly0P + FlyOffsetA);
  1047.             fly4i = *(fly0P + FlyOffsetAIm);
  1048.             fly0r = *(fly1P + FlyOffsetA);
  1049.             fly0i = *(fly1P + FlyOffsetAIm);
  1050.             fly6r = *(fly2P + FlyOffsetA);
  1051.             fly6i = *(fly2P + FlyOffsetAIm);
  1052.             fly3r = *(fly3P + FlyOffsetA);
  1053.             fly3i = *(fly3P + FlyOffsetAIm);
  1054.         
  1055.             fly5r = t0r - fly2r;
  1056.             fly5i = t0i - fly2i;
  1057.             t0r = t0r + fly2r;
  1058.             t0i = t0i + fly2i;
  1059.  
  1060.             fly7r = t1r + fly1i;
  1061.             fly7i = t1i - fly1r;
  1062.             t1r = t1r - fly1i;
  1063.             t1i = t1i + fly1r;
  1064.  
  1065.             fly2r = fly4r - fly6r;
  1066.             fly2i = fly4i - fly6i;
  1067.             fly4r = fly4r + fly6r;
  1068.             fly4i = fly4i + fly6i;
  1069.  
  1070.             fly1r = fly0r + fly3i;
  1071.             fly1i = fly0i - fly3r;
  1072.             fly0r = fly0r - fly3i;
  1073.             fly0i = fly0i + fly3r;
  1074.  
  1075.             *(fly2P) = fly5r;
  1076.             *(fly2P + 1) = fly5i;
  1077.             *(fly0P) = t0r;
  1078.             *(fly0P + 1) = t0i;
  1079.             *(fly3P) = fly7r;
  1080.             *(fly3P + 1) = fly7i;
  1081.             *(fly1P) = t1r;
  1082.             *(fly1P + 1) = t1i;
  1083.             *(fly2P + FlyOffsetA) = fly2r;
  1084.             *(fly2P + FlyOffsetAIm) = fly2i;
  1085.             *(fly0P + FlyOffsetA) = fly4r;
  1086.             *(fly0P + FlyOffsetAIm) = fly4i;
  1087.             *(fly3P + FlyOffsetA) = fly1r;
  1088.             *(fly3P + FlyOffsetAIm) = fly1i;
  1089.             *(fly1P + FlyOffsetA) = fly0r;
  1090.             *(fly1P + FlyOffsetAIm) = fly0i;
  1091.  
  1092.             fly0P = (fly0P + FlyOffsetB);
  1093.             fly1P = (fly1P + FlyOffsetB);
  1094.             fly2P = (fly2P + FlyOffsetB);
  1095.             fly3P = (fly3P + FlyOffsetB);
  1096.         };
  1097.  
  1098.         NsameU4 >>= 1;
  1099.         LoopN >>= 1;
  1100.         NdiffU <<= 1;
  1101.         Flyinc = Flyinc << 1;
  1102.     }
  1103.     else{
  1104.         /* 1 radix 4 stage */
  1105.         IOP = ioptr;
  1106.  
  1107.         U3r =  0.7071067811865475244008443621; /* sqrt(0.5);    */    
  1108.         U3i = U3r;    
  1109.         fly0P = IOP;
  1110.         fly1P = (IOP+Flyinc);
  1111.         fly2P = (fly1P+Flyinc);
  1112.         fly3P = (fly2P+Flyinc);
  1113.         
  1114.             /* Butterflys        */
  1115.             /*
  1116.             t0    -    -    t0    -    -    t0
  1117.             t1    -    -    t1    -    -    t1
  1118.             f2    -  1-    f5    -    -    f5
  1119.             f1    - -i-    f7    -    -    f7
  1120.             f4    -    -    f4    -  1-    f6
  1121.             f0    -    -    f0    -U3    -    f3
  1122.             f6    -  1-    f2    - -i-    f4
  1123.             f3    - -i-    f1    -U3a-    f2
  1124.             */
  1125.  
  1126.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1127.             t0r = *(fly0P);
  1128.             t0i = *(fly0P + 1);
  1129.             t1r = *(fly1P);
  1130.             t1i = *(fly1P + 1);
  1131.             fly2r = *(fly2P);
  1132.             fly2i = *(fly2P + 1);
  1133.             fly1r = *(fly3P);
  1134.             fly1i = *(fly3P + 1);
  1135.             fly4r = *(fly0P + FlyOffsetA);
  1136.             fly4i = *(fly0P + FlyOffsetAIm);
  1137.             fly0r = *(fly1P + FlyOffsetA);
  1138.             fly0i = *(fly1P + FlyOffsetAIm);
  1139.             fly6r = *(fly2P + FlyOffsetA);
  1140.             fly6i = *(fly2P + FlyOffsetAIm);
  1141.             fly3r = *(fly3P + FlyOffsetA);
  1142.             fly3i = *(fly3P + FlyOffsetAIm);
  1143.     
  1144.             fly5r = t0r - fly2r;
  1145.             fly5i = t0i - fly2i;
  1146.             t0r = t0r + fly2r;
  1147.             t0i = t0i + fly2i;
  1148.     
  1149.             fly7r = t1r + fly1i;
  1150.             fly7i = t1i - fly1r;
  1151.             t1r = t1r - fly1i;
  1152.             t1i = t1i + fly1r;
  1153.     
  1154.             fly2r = fly4r - fly6r;
  1155.             fly2i = fly4i - fly6i;
  1156.             fly4r = fly4r + fly6r;
  1157.             fly4i = fly4i + fly6i;
  1158.     
  1159.             fly1r = fly0r + fly3i;
  1160.             fly1i = fly0i - fly3r;
  1161.             fly0r = fly0r - fly3i;
  1162.             fly0i = fly0i + fly3r;
  1163.     
  1164.             fly6r = t0r - fly4r;
  1165.             fly6i = t0i - fly4i;
  1166.             t0r = t0r + fly4r;
  1167.             t0i = t0i + fly4i;
  1168.     
  1169.             fly3r = fly5r + fly2i;
  1170.             fly3i = fly5i - fly2r;
  1171.             fly5r = fly5r - fly2i;
  1172.             fly5i = fly5i + fly2r;
  1173.     
  1174.             fly4r = t1r - U3r * fly0r;
  1175.             fly4r = fly4r + U3i * fly0i;
  1176.             fly4i = t1i - U3i * fly0r;
  1177.             fly4i = fly4i - U3r * fly0i;
  1178.             t1r = scale * t1r - fly4r;
  1179.             t1i = scale * t1i - fly4i;
  1180.     
  1181.             fly2r = fly7r + U3i * fly1r;
  1182.             fly2r = fly2r + U3r * fly1i;
  1183.             fly2i = fly7i - U3r * fly1r;
  1184.             fly2i = fly2i + U3i * fly1i;
  1185.             fly7r = scale * fly7r - fly2r;
  1186.             fly7i = scale * fly7i - fly2i;
  1187.     
  1188.             *(fly0P + FlyOffsetA) = fly6r;
  1189.             *(fly0P + FlyOffsetAIm) = fly6i;
  1190.             *(fly0P) = t0r;
  1191.             *(fly0P + 1) = t0i;
  1192.             *(fly2P + FlyOffsetA) = fly3r;
  1193.             *(fly2P + FlyOffsetAIm) = fly3i;
  1194.             *(fly2P) = fly5r;
  1195.             *(fly2P + 1) = fly5i;
  1196.             *(fly1P + FlyOffsetA) = fly4r;
  1197.             *(fly1P + FlyOffsetAIm) = fly4i;
  1198.             *(fly1P) = t1r;
  1199.             *(fly1P + 1) = t1i;
  1200.             *(fly3P + FlyOffsetA) = fly2r;
  1201.             *(fly3P + FlyOffsetAIm) = fly2i;
  1202.             *(fly3P) = fly7r;
  1203.             *(fly3P + 1) = fly7i;
  1204.         
  1205.             fly0P = (fly0P + FlyOffsetB);
  1206.             fly1P = (fly1P + FlyOffsetB);
  1207.             fly2P = (fly2P + FlyOffsetB);
  1208.             fly3P = (fly3P + FlyOffsetB);
  1209.  
  1210.         };
  1211.  
  1212.         NsameU4 >>= 2;
  1213.         LoopN >>= 2;
  1214.         NdiffU <<= 2;
  1215.         Flyinc = Flyinc << 2;
  1216.     };
  1217. };
  1218.  
  1219. NsameU2 = NsameU4 >> 1;
  1220. NsameU1 = NsameU2 >> 1;
  1221. Flyinc <<= 1;
  1222. FlyOffsetA =  Flyinc << 2;
  1223. FlyOffsetB =  FlyOffsetA << 1;
  1224. FlyOffsetAIm =  FlyOffsetA + 1;
  1225. FlyOffsetBIm =  FlyOffsetB + 1;
  1226. LoopN >>= 1;
  1227.  
  1228. /*   ****** RADIX 8 Stages    */
  1229. for (stage = stage<<1; stage > 0 ; stage--){
  1230.  
  1231.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  1232.  
  1233. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  1234.     if(!(stage&1)){
  1235.         U0rP = &Utbl[0];
  1236.         U0iP = &Utbl[Ntbl[M-2]];
  1237.         U1rP = U0rP;
  1238.         U1iP = U0iP;
  1239.         U2rP = U0rP;
  1240.         U2iP = U0iP;
  1241.         U3offset = (Ntbl[M]) / 8;
  1242.  
  1243.         IOP = ioptr;
  1244.     
  1245.         U0r =  *U0rP,
  1246.         U0i =  *U0iP;
  1247.         U1r =  *U1rP,
  1248.         U1i =  *U1iP;
  1249.         U2r =  *U2rP,
  1250.         U2i =  *U2iP;
  1251.         U3r =  *( U2rP + U3offset);
  1252.         U3i =  *( U2iP - U3offset);
  1253.     }
  1254.     
  1255.     fly0P = IOP;
  1256.     fly1P = (IOP+Flyinc);
  1257.     fly2P = (fly1P+Flyinc);
  1258.     fly3P = (fly2P+Flyinc);
  1259.  
  1260.     for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
  1261.  
  1262.             /* Butterflys        */
  1263.             /*
  1264.             f0    -    -    t0    -    -    t0    -    -    t0
  1265.             f1    -U0    -    t1    -    -    t1    -    -    t1
  1266.             f2    -    -    f2    -U1    -    f5    -    -    f3
  1267.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  1268.             f4    -    -    f4    -    -    f4    -U2    -    f6
  1269.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  1270.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  1271.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  1272.             */
  1273.         
  1274.         fly0r = *(IOP);
  1275.         fly0i = *(IOP+1);
  1276.         fly1r = *(fly1P);
  1277.         fly1i = *(fly1P+1);
  1278.  
  1279.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  1280.     
  1281.             fly2r = *(fly2P);
  1282.             fly2i = *(fly2P + 1);
  1283.             fly3r = *(fly3P);
  1284.             fly3i = *(fly3P + 1);
  1285.             fly4r = *(fly0P + FlyOffsetA);
  1286.             fly4i = *(fly0P + FlyOffsetAIm);
  1287.             fly5r = *(fly1P + FlyOffsetA);
  1288.             fly5i = *(fly1P + FlyOffsetAIm);
  1289.             fly6r = *(fly2P + FlyOffsetA);
  1290.             fly6i = *(fly2P + FlyOffsetAIm);
  1291.             fly7r = *(fly3P + FlyOffsetA);
  1292.             fly7i = *(fly3P + FlyOffsetAIm);
  1293.  
  1294.             t1r = fly0r - U0r * fly1r;
  1295.             t1r = t1r + U0i * fly1i;
  1296.             t1i = fly0i - U0i * fly1r;
  1297.             t1i = t1i - U0r * fly1i;
  1298.             t0r = scale * fly0r - t1r;
  1299.             t0i = scale * fly0i - t1i;
  1300.     
  1301.             fly1r = fly2r - U0r * fly3r;
  1302.             fly1r = fly1r + U0i * fly3i;
  1303.             fly1i = fly2i - U0i * fly3r;
  1304.             fly1i = fly1i - U0r * fly3i;
  1305.             fly2r = scale * fly2r - fly1r;
  1306.             fly2i = scale * fly2i - fly1i;
  1307.     
  1308.             fly0r = fly4r - U0r * fly5r;
  1309.             fly0r = fly0r + U0i * fly5i;
  1310.             fly0i = fly4i - U0i * fly5r;
  1311.             fly0i = fly0i - U0r * fly5i;
  1312.             fly4r = scale * fly4r - fly0r;
  1313.             fly4i = scale * fly4i - fly0i;
  1314.     
  1315.             fly3r = fly6r - U0r * fly7r;
  1316.             fly3r = fly3r + U0i * fly7i;
  1317.             fly3i = fly6i - U0i * fly7r;
  1318.             fly3i = fly3i - U0r * fly7i;
  1319.             fly6r = scale * fly6r - fly3r;
  1320.             fly6i = scale * fly6i - fly3i;
  1321.     
  1322.  
  1323.             fly5r = t0r - U1r * fly2r;
  1324.             fly5r = fly5r + U1i * fly2i;
  1325.             fly5i = t0i - U1i * fly2r;
  1326.             fly5i = fly5i - U1r * fly2i;
  1327.             t0r = scale * t0r - fly5r;
  1328.             t0i = scale * t0i - fly5i;
  1329.  
  1330.             fly7r = t1r + U1i * fly1r;
  1331.             fly7r = fly7r + U1r * fly1i;
  1332.             fly7i = t1i - U1r * fly1r;
  1333.             fly7i = fly7i + U1i * fly1i;
  1334.             t1r = scale * t1r - fly7r;
  1335.             t1i = scale * t1i - fly7i;
  1336.  
  1337.             fly2r = fly4r - U1r * fly6r;
  1338.             fly2r = fly2r + U1i * fly6i;
  1339.             fly2i = fly4i - U1i * fly6r;
  1340.             fly2i = fly2i - U1r * fly6i;
  1341.             fly4r = scale * fly4r - fly2r;
  1342.             fly4i = scale * fly4i - fly2i;
  1343.  
  1344.             fly1r = fly0r + U1i * fly3r;
  1345.             fly1r = fly1r + U1r * fly3i;
  1346.             fly1i = fly0i - U1r * fly3r;
  1347.             fly1i = fly1i + U1i * fly3i;
  1348.             fly0r = scale * fly0r - fly1r;
  1349.             fly0i = scale * fly0i - fly1i;
  1350.  
  1351.             fly6r = t0r - U2r * fly4r;
  1352.             fly6r = fly6r + U2i * fly4i;
  1353.             fly6i = t0i - U2i * fly4r;
  1354.             fly6i = fly6i - U2r * fly4i;
  1355.             t0r = scale * t0r - fly6r;
  1356.             t0i = scale * t0i - fly6i;
  1357.  
  1358.             fly3r = fly5r - U2i * fly2r;
  1359.             fly3r = fly3r - U2r * fly2i;
  1360.             fly3i = fly5i + U2r * fly2r;
  1361.             fly3i = fly3i - U2i * fly2i;
  1362.             fly2r = scale * fly5r - fly3r;
  1363.             fly2i = scale * fly5i - fly3i;
  1364.  
  1365.             fly4r = t1r - U3r * fly0r;
  1366.             fly4r = fly4r + U3i * fly0i;
  1367.             fly4i = t1i - U3i * fly0r;
  1368.             fly4i = fly4i - U3r * fly0i;
  1369.             t1r = scale * t1r - fly4r;
  1370.             t1i = scale * t1i - fly4i;
  1371.  
  1372.             fly5r = fly7r + U3i * fly1r;
  1373.             fly5r = fly5r + U3r * fly1i;
  1374.             fly5i = fly7i - U3r * fly1r;
  1375.             fly5i = fly5i + U3i * fly1i;
  1376.             fly7r = scale * fly7r - fly5r;
  1377.             fly7i = scale * fly7i - fly5i;
  1378.  
  1379.             *(fly0P + FlyOffsetA) = fly6r;
  1380.             *(fly0P + FlyOffsetAIm) = fly6i;
  1381.             *(fly0P) = t0r;
  1382.             *(fly0P + 1) = t0i;
  1383.             *(fly2P) = fly3r;
  1384.             *(fly2P + 1) = fly3i;
  1385.             *(fly2P + FlyOffsetA) = fly2r;
  1386.             *(fly2P + FlyOffsetAIm) = fly2i;
  1387.  
  1388.             fly0r = *(fly0P + FlyOffsetB);
  1389.             fly0i = *(fly0P + FlyOffsetBIm);
  1390.  
  1391.             *(fly1P + FlyOffsetA) = fly4r;
  1392.             *(fly1P + FlyOffsetAIm) = fly4i;
  1393.             *(fly1P) = t1r;
  1394.             *(fly1P + 1) = t1i;
  1395.  
  1396.             fly1r = *(fly1P + FlyOffsetB);
  1397.             fly1i = *(fly1P + FlyOffsetBIm);
  1398.  
  1399.             *(fly3P + FlyOffsetA) = fly5r;
  1400.             *(fly3P + FlyOffsetAIm) = fly5i;
  1401.             *(fly3P) = fly7r;
  1402.             *(fly3P + 1) = fly7i;
  1403.  
  1404.             fly0P = (fly0P + FlyOffsetB);
  1405.             fly1P = (fly1P + FlyOffsetB);
  1406.             fly2P = (fly2P + FlyOffsetB);
  1407.             fly3P = (fly3P + FlyOffsetB);
  1408.  
  1409.         };
  1410.         fly2r = *(fly2P);
  1411.         fly2i = *(fly2P + 1);
  1412.         fly3r = *(fly3P);
  1413.         fly3i = *(fly3P + 1);
  1414.         fly4r = *(fly0P + FlyOffsetA);
  1415.         fly4i = *(fly0P + FlyOffsetAIm);
  1416.         fly5r = *(fly1P + FlyOffsetA);
  1417.         fly5i = *(fly1P + FlyOffsetAIm);
  1418.         fly6r = *(fly2P + FlyOffsetA);
  1419.         fly6i = *(fly2P + FlyOffsetAIm);
  1420.         fly7r = *(fly3P + FlyOffsetA);
  1421.         fly7i = *(fly3P + FlyOffsetAIm);
  1422.  
  1423.         t1r = fly0r - U0r * fly1r;
  1424.         t1r = t1r + U0i * fly1i;
  1425.         t1i = fly0i - U0i * fly1r;
  1426.         t1i = t1i - U0r * fly1i;
  1427.         t0r = scale * fly0r - t1r;
  1428.         t0i = scale * fly0i - t1i;
  1429.  
  1430.         fly1r = fly2r - U0r * fly3r;
  1431.         fly1r = fly1r + U0i * fly3i;
  1432.         fly1i = fly2i - U0i * fly3r;
  1433.         fly1i = fly1i - U0r * fly3i;
  1434.         fly2r = scale * fly2r - fly1r;
  1435.         fly2i = scale * fly2i - fly1i;
  1436.  
  1437.         fly0r = fly4r - U0r * fly5r;
  1438.         fly0r = fly0r + U0i * fly5i;
  1439.         fly0i = fly4i - U0i * fly5r;
  1440.         fly0i = fly0i - U0r * fly5i;
  1441.         fly4r = scale * fly4r - fly0r;
  1442.         fly4i = scale * fly4i - fly0i;
  1443.  
  1444.         fly3r = fly6r - U0r * fly7r;
  1445.         fly3r = fly3r + U0i * fly7i;
  1446.         fly3i = fly6i - U0i * fly7r;
  1447.         fly3i = fly3i - U0r * fly7i;
  1448.         fly6r = scale * fly6r - fly3r;
  1449.         fly6i = scale * fly6i - fly3i;
  1450.  
  1451.         fly5r = t0r - U1r * fly2r;
  1452.         fly5r = fly5r + U1i * fly2i;
  1453.         fly5i = t0i - U1i * fly2r;
  1454.         fly5i = fly5i - U1r * fly2i;
  1455.         t0r = scale * t0r - fly5r;
  1456.         t0i = scale * t0i - fly5i;
  1457.  
  1458.         fly7r = t1r + U1i * fly1r;
  1459.         fly7r = fly7r + U1r * fly1i;
  1460.         fly7i = t1i - U1r * fly1r;
  1461.         fly7i = fly7i + U1i * fly1i;
  1462.         t1r = scale * t1r - fly7r;
  1463.         t1i = scale * t1i - fly7i;
  1464.  
  1465.         fly2r = fly4r - U1r * fly6r;
  1466.         fly2r = fly2r + U1i * fly6i;
  1467.         fly2i = fly4i - U1i * fly6r;
  1468.         fly2i = fly2i - U1r * fly6i;
  1469.         fly4r = scale * fly4r - fly2r;
  1470.         fly4i = scale * fly4i - fly2i;
  1471.  
  1472.         fly1r = fly0r + U1i * fly3r;
  1473.         fly1r = fly1r + U1r * fly3i;
  1474.         fly1i = fly0i - U1r * fly3r;
  1475.         fly1i = fly1i + U1i * fly3i;
  1476.         fly0r = scale * fly0r - fly1r;
  1477.         fly0i = scale * fly0i - fly1i;
  1478.  
  1479.         U0i = *(U0iP = (U0iP - NsameU4));
  1480.         U0r = *(U0rP = (U0rP + NsameU4));        
  1481.         U1r = *(U1rP = (U1rP + NsameU2));
  1482.         U1i = *(U1iP = (U1iP - NsameU2));
  1483.         if(stage&1)
  1484.             U0r = - U0r;
  1485.  
  1486.         fly6r = t0r - U2r * fly4r;
  1487.         fly6r = fly6r + U2i * fly4i;
  1488.         fly6i = t0i - U2i * fly4r;
  1489.         fly6i = fly6i - U2r * fly4i;
  1490.         t0r = scale * t0r - fly6r;
  1491.         t0i = scale * t0i - fly6i;
  1492.  
  1493.         fly3r = fly5r - U2i * fly2r;
  1494.         fly3r = fly3r - U2r * fly2i;
  1495.         fly3i = fly5i + U2r * fly2r;
  1496.         fly3i = fly3i - U2i * fly2i;
  1497.         fly2r = scale * fly5r - fly3r;
  1498.         fly2i = scale * fly5i - fly3i;
  1499.  
  1500.         fly4r = t1r - U3r * fly0r;
  1501.         fly4r = fly4r + U3i * fly0i;
  1502.         fly4i = t1i - U3i * fly0r;
  1503.         fly4i = fly4i - U3r * fly0i;
  1504.         t1r = scale * t1r - fly4r;
  1505.         t1i = scale * t1i - fly4i;
  1506.  
  1507.         fly5r = fly7r + U3i * fly1r;
  1508.         fly5r = fly5r + U3r * fly1i;
  1509.         fly5i = fly7i - U3r * fly1r;
  1510.         fly5i = fly5i + U3i * fly1i;
  1511.         fly7r = scale * fly7r - fly5r;
  1512.         fly7i = scale * fly7i - fly5i;
  1513.  
  1514.         *(fly0P + FlyOffsetA) = fly6r;
  1515.         *(fly0P + FlyOffsetAIm) = fly6i;
  1516.         *(fly0P) = t0r;
  1517.         *(fly0P + 1) = t0i;
  1518.  
  1519.         U2r = *(U2rP = (U2rP + NsameU1));
  1520.         U2i = *(U2iP = (U2iP - NsameU1));
  1521.  
  1522.         *(fly2P) = fly3r;
  1523.         *(fly2P + 1) = fly3i;
  1524.         *(fly2P + FlyOffsetA) = fly2r;
  1525.         *(fly2P + FlyOffsetAIm) = fly2i;
  1526.         *(fly1P + FlyOffsetA) = fly4r;
  1527.         *(fly1P + FlyOffsetAIm) = fly4i;
  1528.         *(fly1P) = t1r;
  1529.         *(fly1P + 1) = t1i;
  1530.  
  1531.         U3r = *( U2rP + U3offset);
  1532.         U3i = *( U2iP - U3offset);
  1533.  
  1534.         *(fly3P + FlyOffsetA) = fly5r;
  1535.         *(fly3P + FlyOffsetAIm) = fly5i;
  1536.         *(fly3P) = fly7r;
  1537.         *(fly3P + 1) = fly7i;
  1538.  
  1539.         IOP = IOP + 2;
  1540.         fly0P = IOP;
  1541.         fly1P = (IOP+Flyinc);
  1542.         fly2P = (fly1P+Flyinc);
  1543.         fly3P = (fly2P+Flyinc);
  1544.     };
  1545.     NsameU4 = -NsameU4;
  1546.  
  1547.     if(stage&1){
  1548.         LoopN >>= 3;
  1549.         NsameU1 >>= 3;
  1550.         NsameU2 >>= 3;
  1551.         NsameU4 >>= 3;
  1552.         NdiffU <<= 3;
  1553.         Flyinc = Flyinc << 3;
  1554.         FlyOffsetA <<= 3;
  1555.         FlyOffsetB <<= 3;
  1556.         FlyOffsetAIm =  FlyOffsetA + 1;
  1557.         FlyOffsetBIm =  FlyOffsetB + 1;
  1558.     }
  1559. }
  1560. ioptr += 2*Ntbl[M];
  1561. }
  1562. }
  1563.  
  1564. /* rffts    */
  1565. /* compute multiple real input FFTs on 'Rows' consecutively stored arrays    */
  1566. /* ioptr = pointer to the data in and out    */
  1567. /* M = log2 of fft size    */
  1568. /* Rows = number of FFTs to compute    */
  1569. /* Utbl = Pointer to cosine table    */
  1570.  
  1571. void rffts(float *ioptr, long M, long Rows, float *Utbl){
  1572. /* Compute real fft on the rows of the input array    */
  1573. /* INPUTS */
  1574. /* M = log2 of fft size    */
  1575. /* *ioptr = input data array    */
  1576. /* *Utbl = cosine table    */
  1577. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  1578. /* OUTPUTS */
  1579. /* *ioptr = output data array    */
  1580.  
  1581. long     Flyinc;
  1582. long    FlyOffsetA;
  1583. long    FlyOffsetAIm;
  1584. long    FlyOffsetB;
  1585. long    FlyOffsetBIm;
  1586. long     NsameU1;
  1587. long     NsameU2;
  1588. long     NsameU4;
  1589. long     diffUcnt;
  1590. long     LoopCnt;
  1591. float    scale;
  1592. float    fly0r;
  1593. float    fly0i;
  1594. float    fly1r;
  1595. float    fly1i;
  1596. float    fly2r;
  1597. float    fly2i;
  1598. float    fly3r;
  1599. float    fly3i;
  1600. float    fly4r;
  1601. float    fly4i;
  1602. float    fly5r;
  1603. float    fly5i;
  1604. float    fly6r;
  1605. float    fly6i;
  1606. float    fly7r;
  1607. float    fly7i;
  1608. float    U0r;
  1609. float    U0i;
  1610. float    U1r;
  1611. float    U1i;
  1612. float    U2r;
  1613. float    U2i;
  1614. float    U3r;
  1615. float    U3i;
  1616. float    t0r;
  1617. float    t0i;
  1618. float    t1r;
  1619. float    t1i;
  1620.  
  1621. float    *fly0P;
  1622. float    *fly1P;
  1623. float    *fly2P;
  1624. float    *fly3P;
  1625.  
  1626. float    *U0rP;
  1627. float    *U0iP;
  1628. float    *U1rP;
  1629. float    *U1iP;
  1630. float    *U2rP;
  1631. float    *U2iP;
  1632. float    *IOP;
  1633. long    U3offset;
  1634.  
  1635. long     stage;
  1636. long     NdiffU;
  1637. long     LoopN;
  1638.  
  1639. const long BRshift = MAXMROOT - ((M-1)>>1);        /* for RFFT */
  1640. const long Nrems2 = Ntbl[(M-1)-((M-1)>>1)+1];        /* for RFFT */
  1641. const long Nroot_1_ColInc = (Ntbl[(M-1)-1]-Ntbl[(M-1)-((M-1)>>1)])*2; /* for RFFT */
  1642.  
  1643. for (;Rows>0;Rows--){
  1644.  
  1645. M=M-1;        /* for RFFT */
  1646.  
  1647. FlyOffsetA = Ntbl[M] * 2/2;
  1648. FlyOffsetAIm = FlyOffsetA + 1;
  1649. FlyOffsetB = FlyOffsetA + 2;
  1650. FlyOffsetBIm = FlyOffsetB + 1;
  1651.  
  1652. /* BitrevR2 ** bit reverse shuffle and first radix 2 stage ******/
  1653.  
  1654. scale = 0.5;
  1655. for (stage = 0; stage < Ntbl[M-(M>>1)]*2; stage += Ntbl[M>>1]*2){
  1656.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  1657.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  1658.         fly0P = ioptr + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(2*2) + stage;
  1659.         IOP = ioptr + (LoopN<<(M+1)/2) * 2 + stage;
  1660.         fly1P = IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2);
  1661.         fly0r = *fly0P;
  1662.         fly0i = *(fly0P+1);
  1663.         fly1r = *(fly0P+FlyOffsetA);
  1664.         fly1i = *(fly0P+FlyOffsetAIm);
  1665.         for (; LoopCnt > LoopN;){
  1666.             fly2r = *(fly0P+2);
  1667.             fly2i = *(fly0P+(2+1));
  1668.             fly3r = *(fly0P+FlyOffsetB);
  1669.             fly3i = *(fly0P+FlyOffsetBIm);
  1670.             fly4r = *fly1P;
  1671.             fly4i = *(fly1P+1);
  1672.             fly5r = *(fly1P+FlyOffsetA);
  1673.             fly5i = *(fly1P+FlyOffsetAIm);
  1674.             fly6r = *(fly1P+2);
  1675.             fly6i = *(fly1P+(2+1));
  1676.             fly7r = *(fly1P+FlyOffsetB);
  1677.             fly7i = *(fly1P+FlyOffsetBIm);
  1678.             
  1679.             t0r    = fly0r + fly1r;
  1680.             t0i    = fly0i + fly1i;
  1681.             fly1r = fly0r - fly1r;
  1682.             fly1i = fly0i - fly1i;
  1683.             t1r = fly2r + fly3r;
  1684.             t1i = fly2i + fly3i;
  1685.             fly3r = fly2r - fly3r;
  1686.             fly3i = fly2i - fly3i;
  1687.             fly0r = fly4r + fly5r;
  1688.             fly0i = fly4i + fly5i;
  1689.             fly5r = fly4r - fly5r;
  1690.             fly5i = fly4i - fly5i;
  1691.             fly2r = fly6r + fly7r;
  1692.             fly2i = fly6i + fly7i;
  1693.             fly7r = fly6r - fly7r;
  1694.             fly7i = fly6i - fly7i;
  1695.  
  1696.             *fly1P = scale*t0r;
  1697.             *(fly1P+1) = scale*t0i;
  1698.             *(fly1P+2) = scale*fly1r;
  1699.             *(fly1P+(2+1)) = scale*fly1i;
  1700.             *(fly1P+FlyOffsetA) = scale*t1r;
  1701.             *(fly1P+FlyOffsetAIm) = scale*t1i;
  1702.             *(fly1P+FlyOffsetB) = scale*fly3r;
  1703.             *(fly1P+FlyOffsetBIm) = scale*fly3i;
  1704.             *fly0P = scale*fly0r;
  1705.             *(fly0P+1) = scale*fly0i;
  1706.             *(fly0P+2) = scale*fly5r;
  1707.             *(fly0P+(2+1)) = scale*fly5i;
  1708.             *(fly0P+FlyOffsetA) = scale*fly2r;
  1709.             *(fly0P+FlyOffsetAIm) = scale*fly2i;
  1710.             *(fly0P+FlyOffsetB) = scale*fly7r;
  1711.             *(fly0P+FlyOffsetBIm) = scale*fly7i;
  1712.  
  1713.             fly0P -= Nrems2;
  1714.             fly0r = *fly0P;
  1715.             fly0i = *(fly0P+1);
  1716.             fly1r = *(fly0P+FlyOffsetA);
  1717.             fly1i = *(fly0P+FlyOffsetAIm);
  1718.             LoopCnt -= 1;
  1719.             fly1P = (IOP + ((int)BRcnt[LoopCnt] >> BRshift)*(2*2));
  1720.         };
  1721.         fly2r = *(fly0P+2);
  1722.         fly2i = *(fly0P+(2+1));
  1723.         fly3r = *(fly0P+FlyOffsetB);
  1724.         fly3i = *(fly0P+FlyOffsetBIm);
  1725.  
  1726.         t0r   = fly0r + fly1r;
  1727.         t0i   = fly0i + fly1i;
  1728.         fly1r = fly0r - fly1r;
  1729.         fly1i = fly0i - fly1i;
  1730.         t1r   = fly2r + fly3r;
  1731.         t1i   = fly2i + fly3i;
  1732.         fly3r = fly2r - fly3r;
  1733.         fly3i = fly2i - fly3i;
  1734.  
  1735.         *fly0P = scale*t0r;
  1736.         *(fly0P+1) = scale*t0i;
  1737.         *(fly0P+2) = scale*fly1r;
  1738.         *(fly0P+(2+1)) = scale*fly1i;
  1739.         *(fly0P+FlyOffsetA) = scale*t1r;
  1740.         *(fly0P+FlyOffsetAIm) = scale*t1i;
  1741.         *(fly0P+FlyOffsetB) = scale*fly3r;
  1742.         *(fly0P+FlyOffsetBIm) = scale*fly3i;
  1743.  
  1744.     };
  1745. };
  1746.  
  1747.  
  1748.  
  1749. /**** FFTC  **************/
  1750.  
  1751. scale = 2.0;
  1752.  
  1753. NdiffU = 2;
  1754. Flyinc =  (NdiffU);
  1755.  
  1756. NsameU4 = Ntbl[M+1]/4;    /* for RFFT */
  1757. LoopN = Ntbl[M-3];
  1758.  
  1759. stage = ((M-1)/3);
  1760. if ( (M-1-(stage * 3)) != 0 ){
  1761.     FlyOffsetA =  Flyinc << 2;
  1762.     FlyOffsetB =  FlyOffsetA << 1;
  1763.     FlyOffsetAIm =  FlyOffsetA + 1;
  1764.     FlyOffsetBIm =  FlyOffsetB + 1;
  1765.     if ( (M-1-(stage * 3)) == 1 ){
  1766.         /* 1 radix 2 stage */
  1767.  
  1768.         IOP = ioptr;
  1769.         fly0P = IOP;
  1770.         fly1P = (IOP+Flyinc);
  1771.         fly2P = (fly1P+Flyinc);
  1772.         fly3P = (fly2P+Flyinc);
  1773.         
  1774.             /* Butterflys        */
  1775.             /*
  1776.             t0    -    -    t0
  1777.             t1    -    -    t1
  1778.             f2    -  1-    f5
  1779.             f1    - -i-    f7
  1780.             f4    -    -    f4
  1781.             f0    -    -    f0
  1782.             f6    -  1-    f2
  1783.             f3    - -i-    f1
  1784.             */
  1785.  
  1786.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1787.             t0r = *(fly0P);
  1788.             t0i = *(fly0P + 1);
  1789.             t1r = *(fly1P);
  1790.             t1i = *(fly1P + 1);
  1791.             fly2r = *(fly2P);
  1792.             fly2i = *(fly2P + 1);
  1793.             fly1r = *(fly3P);
  1794.             fly1i = *(fly3P + 1);
  1795.             fly4r = *(fly0P + FlyOffsetA);
  1796.             fly4i = *(fly0P + FlyOffsetAIm);
  1797.             fly0r = *(fly1P + FlyOffsetA);
  1798.             fly0i = *(fly1P + FlyOffsetAIm);
  1799.             fly6r = *(fly2P + FlyOffsetA);
  1800.             fly6i = *(fly2P + FlyOffsetAIm);
  1801.             fly3r = *(fly3P + FlyOffsetA);
  1802.             fly3i = *(fly3P + FlyOffsetAIm);
  1803.         
  1804.             fly5r = t0r - fly2r;
  1805.             fly5i = t0i - fly2i;
  1806.             t0r = t0r + fly2r;
  1807.             t0i = t0i + fly2i;
  1808.  
  1809.             fly7r = t1r - fly1i;
  1810.             fly7i = t1i + fly1r;
  1811.             t1r = t1r + fly1i;
  1812.             t1i = t1i - fly1r;
  1813.  
  1814.             fly2r = fly4r - fly6r;
  1815.             fly2i = fly4i - fly6i;
  1816.             fly4r = fly4r + fly6r;
  1817.             fly4i = fly4i + fly6i;
  1818.  
  1819.             fly1r = fly0r - fly3i;
  1820.             fly1i = fly0i + fly3r;
  1821.             fly0r = fly0r + fly3i;
  1822.             fly0i = fly0i - fly3r;
  1823.  
  1824.             *(fly2P) = fly5r;
  1825.             *(fly2P + 1) = fly5i;
  1826.             *(fly0P) = t0r;
  1827.             *(fly0P + 1) = t0i;
  1828.             *(fly3P) = fly7r;
  1829.             *(fly3P + 1) = fly7i;
  1830.             *(fly1P) = t1r;
  1831.             *(fly1P + 1) = t1i;
  1832.             *(fly2P + FlyOffsetA) = fly2r;
  1833.             *(fly2P + FlyOffsetAIm) = fly2i;
  1834.             *(fly0P + FlyOffsetA) = fly4r;
  1835.             *(fly0P + FlyOffsetAIm) = fly4i;
  1836.             *(fly3P + FlyOffsetA) = fly1r;
  1837.             *(fly3P + FlyOffsetAIm) = fly1i;
  1838.             *(fly1P + FlyOffsetA) = fly0r;
  1839.             *(fly1P + FlyOffsetAIm) = fly0i;
  1840.  
  1841.             fly0P = (fly0P + FlyOffsetB);
  1842.             fly1P = (fly1P + FlyOffsetB);
  1843.             fly2P = (fly2P + FlyOffsetB);
  1844.             fly3P = (fly3P + FlyOffsetB);
  1845.         };
  1846.  
  1847.         NsameU4 >>= 1;
  1848.         LoopN >>= 1;
  1849.         NdiffU <<= 1;
  1850.         Flyinc = Flyinc << 1;
  1851.     }
  1852.     else{
  1853.         /* 1 radix 4 stage */
  1854.         IOP = ioptr;
  1855.  
  1856.         U3r =  0.7071067811865475244008443621; /* sqrt(0.5);    */    
  1857.         U3i = U3r;    
  1858.         fly0P = IOP;
  1859.         fly1P = (IOP+Flyinc);
  1860.         fly2P = (fly1P+Flyinc);
  1861.         fly3P = (fly2P+Flyinc);
  1862.         
  1863.             /* Butterflys        */
  1864.             /*
  1865.             t0    -    -    t0    -    -    t0
  1866.             t1    -    -    t1    -    -    t1
  1867.             f2    -  1-    f5    -    -    f5
  1868.             f1    - -i-    f7    -    -    f7
  1869.             f4    -    -    f4    -  1-    f6
  1870.             f0    -    -    f0    -U3    -    f3
  1871.             f6    -  1-    f2    - -i-    f4
  1872.             f3    - -i-    f1    -U3a-    f2
  1873.             */
  1874.  
  1875.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1876.             t0r = *(fly0P);
  1877.             t0i = *(fly0P + 1);
  1878.             t1r = *(fly1P);
  1879.             t1i = *(fly1P + 1);
  1880.             fly2r = *(fly2P);
  1881.             fly2i = *(fly2P + 1);
  1882.             fly1r = *(fly3P);
  1883.             fly1i = *(fly3P + 1);
  1884.             fly4r = *(fly0P + FlyOffsetA);
  1885.             fly4i = *(fly0P + FlyOffsetAIm);
  1886.             fly0r = *(fly1P + FlyOffsetA);
  1887.             fly0i = *(fly1P + FlyOffsetAIm);
  1888.             fly6r = *(fly2P + FlyOffsetA);
  1889.             fly6i = *(fly2P + FlyOffsetAIm);
  1890.             fly3r = *(fly3P + FlyOffsetA);
  1891.             fly3i = *(fly3P + FlyOffsetAIm);
  1892.     
  1893.             fly5r = t0r - fly2r;
  1894.             fly5i = t0i - fly2i;
  1895.             t0r = t0r + fly2r;
  1896.             t0i = t0i + fly2i;
  1897.     
  1898.             fly7r = t1r - fly1i;
  1899.             fly7i = t1i + fly1r;
  1900.             t1r = t1r + fly1i;
  1901.             t1i = t1i - fly1r;
  1902.     
  1903.             fly2r = fly4r - fly6r;
  1904.             fly2i = fly4i - fly6i;
  1905.             fly4r = fly4r + fly6r;
  1906.             fly4i = fly4i + fly6i;
  1907.     
  1908.             fly1r = fly0r - fly3i;
  1909.             fly1i = fly0i + fly3r;
  1910.             fly0r = fly0r + fly3i;
  1911.             fly0i = fly0i - fly3r;
  1912.     
  1913.             fly6r = t0r - fly4r;
  1914.             fly6i = t0i - fly4i;
  1915.             t0r = t0r + fly4r;
  1916.             t0i = t0i + fly4i;
  1917.     
  1918.             fly3r = fly5r - fly2i;
  1919.             fly3i = fly5i + fly2r;
  1920.             fly5r = fly5r + fly2i;
  1921.             fly5i = fly5i - fly2r;
  1922.     
  1923.             fly4r = t1r - U3r * fly0r;
  1924.             fly4r = fly4r - U3i * fly0i;
  1925.             fly4i = t1i + U3i * fly0r;
  1926.             fly4i = fly4i - U3r * fly0i;
  1927.             t1r = scale * t1r - fly4r;
  1928.             t1i = scale * t1i - fly4i;
  1929.     
  1930.             fly2r = fly7r + U3i * fly1r;
  1931.             fly2r = fly2r - U3r * fly1i;
  1932.             fly2i = fly7i + U3r * fly1r;
  1933.             fly2i = fly2i + U3i * fly1i;
  1934.             fly7r = scale * fly7r - fly2r;
  1935.             fly7i = scale * fly7i - fly2i;
  1936.     
  1937.             *(fly0P + FlyOffsetA) = fly6r;
  1938.             *(fly0P + FlyOffsetAIm) = fly6i;
  1939.             *(fly0P) = t0r;
  1940.             *(fly0P + 1) = t0i;
  1941.             *(fly2P + FlyOffsetA) = fly3r;
  1942.             *(fly2P + FlyOffsetAIm) = fly3i;
  1943.             *(fly2P) = fly5r;
  1944.             *(fly2P + 1) = fly5i;
  1945.             *(fly1P + FlyOffsetA) = fly4r;
  1946.             *(fly1P + FlyOffsetAIm) = fly4i;
  1947.             *(fly1P) = t1r;
  1948.             *(fly1P + 1) = t1i;
  1949.             *(fly3P + FlyOffsetA) = fly2r;
  1950.             *(fly3P + FlyOffsetAIm) = fly2i;
  1951.             *(fly3P) = fly7r;
  1952.             *(fly3P + 1) = fly7i;
  1953.         
  1954.             fly0P = (fly0P + FlyOffsetB);
  1955.             fly1P = (fly1P + FlyOffsetB);
  1956.             fly2P = (fly2P + FlyOffsetB);
  1957.             fly3P = (fly3P + FlyOffsetB);
  1958.  
  1959.         };
  1960.         
  1961.         NsameU4 >>= 2;
  1962.         LoopN >>= 2;
  1963.         NdiffU <<= 2;
  1964.         Flyinc = Flyinc << 2;
  1965.     };
  1966. };
  1967.  
  1968. NsameU2 = NsameU4 >> 1;
  1969. NsameU1 = NsameU2 >> 1;
  1970. Flyinc <<= 1;
  1971. FlyOffsetA =  Flyinc << 2;
  1972. FlyOffsetB =  FlyOffsetA << 1;
  1973. FlyOffsetAIm =  FlyOffsetA + 1;
  1974. FlyOffsetBIm =  FlyOffsetB + 1;
  1975. LoopN >>= 1;
  1976.  
  1977. /*   ****** RADIX 8 Stages    */
  1978. for (stage = stage<<1; stage > 0 ; stage--){
  1979.  
  1980.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  1981.  
  1982. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  1983.     if(!(stage&1)){
  1984.         U0rP = &Utbl[0];
  1985.         U0iP = &Utbl[Ntbl[M-1]];    /* for RFFT */
  1986.         U1rP = U0rP;
  1987.         U1iP = U0iP;
  1988.         U2rP = U0rP;
  1989.         U2iP = U0iP;
  1990.         U3offset = (Ntbl[M+1]) / 8;    /* for RFFT */
  1991.  
  1992.         IOP = ioptr;
  1993.     
  1994.         U0r =  *U0rP,
  1995.         U0i =  *U0iP;
  1996.         U1r =  *U1rP,
  1997.         U1i =  *U1iP;
  1998.         U2r =  *U2rP,
  1999.         U2i =  *U2iP;
  2000.         U3r =  *( U2rP + U3offset);
  2001.         U3i =  *( U2iP - U3offset);
  2002.     }
  2003.     
  2004.     fly0P = IOP;
  2005.     fly1P = (IOP+Flyinc);
  2006.     fly2P = (fly1P+Flyinc);
  2007.     fly3P = (fly2P+Flyinc);
  2008.  
  2009.     for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
  2010.  
  2011.             /* Butterflys        */
  2012.             /*
  2013.             f0    -    -    t0    -    -    t0    -    -    t0
  2014.             f1    -U0    -    t1    -    -    t1    -    -    t1
  2015.             f2    -    -    f2    -U1    -    f5    -    -    f3
  2016.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  2017.             f4    -    -    f4    -    -    f4    -U2    -    f6
  2018.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  2019.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  2020.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  2021.             */
  2022.         
  2023.         fly0r = *(IOP);
  2024.         fly0i = *(IOP+1);
  2025.         fly1r = *(fly1P);
  2026.         fly1i = *(fly1P+1);
  2027.  
  2028.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  2029.     
  2030.             fly2r = *(fly2P);
  2031.             fly2i = *(fly2P + 1);
  2032.             fly3r = *(fly3P);
  2033.             fly3i = *(fly3P + 1);
  2034.             fly4r = *(fly0P + FlyOffsetA);
  2035.             fly4i = *(fly0P + FlyOffsetAIm);
  2036.             fly5r = *(fly1P + FlyOffsetA);
  2037.             fly5i = *(fly1P + FlyOffsetAIm);
  2038.             fly6r = *(fly2P + FlyOffsetA);
  2039.             fly6i = *(fly2P + FlyOffsetAIm);
  2040.             fly7r = *(fly3P + FlyOffsetA);
  2041.             fly7i = *(fly3P + FlyOffsetAIm);
  2042.  
  2043.             t1r = fly0r - U0r * fly1r;
  2044.             t1r = t1r - U0i * fly1i;
  2045.             t1i = fly0i + U0i * fly1r;
  2046.             t1i = t1i - U0r * fly1i;
  2047.             t0r = scale * fly0r - t1r;
  2048.             t0i = scale * fly0i - t1i;
  2049.     
  2050.             fly1r = fly2r - U0r * fly3r;
  2051.             fly1r = fly1r - U0i * fly3i;
  2052.             fly1i = fly2i + U0i * fly3r;
  2053.             fly1i = fly1i - U0r * fly3i;
  2054.             fly2r = scale * fly2r - fly1r;
  2055.             fly2i = scale * fly2i - fly1i;
  2056.     
  2057.             fly0r = fly4r - U0r * fly5r;
  2058.             fly0r = fly0r - U0i * fly5i;
  2059.             fly0i = fly4i + U0i * fly5r;
  2060.             fly0i = fly0i - U0r * fly5i;
  2061.             fly4r = scale * fly4r - fly0r;
  2062.             fly4i = scale * fly4i - fly0i;
  2063.     
  2064.             fly3r = fly6r - U0r * fly7r;
  2065.             fly3r = fly3r - U0i * fly7i;
  2066.             fly3i = fly6i + U0i * fly7r;
  2067.             fly3i = fly3i - U0r * fly7i;
  2068.             fly6r = scale * fly6r - fly3r;
  2069.             fly6i = scale * fly6i - fly3i;
  2070.     
  2071.  
  2072.             fly5r = t0r - U1r * fly2r;
  2073.             fly5r = fly5r - U1i * fly2i;
  2074.             fly5i = t0i + U1i * fly2r;
  2075.             fly5i = fly5i - U1r * fly2i;
  2076.             t0r = scale * t0r - fly5r;
  2077.             t0i = scale * t0i - fly5i;
  2078.  
  2079.             fly7r = t1r + U1i * fly1r;
  2080.             fly7r = fly7r - U1r * fly1i;
  2081.             fly7i = t1i + U1r * fly1r;
  2082.             fly7i = fly7i + U1i * fly1i;
  2083.             t1r = scale * t1r - fly7r;
  2084.             t1i = scale * t1i - fly7i;
  2085.  
  2086.             fly2r = fly4r - U1r * fly6r;
  2087.             fly2r = fly2r - U1i * fly6i;
  2088.             fly2i = fly4i + U1i * fly6r;
  2089.             fly2i = fly2i - U1r * fly6i;
  2090.             fly4r = scale * fly4r - fly2r;
  2091.             fly4i = scale * fly4i - fly2i;
  2092.  
  2093.             fly1r = fly0r + U1i * fly3r;
  2094.             fly1r = fly1r - U1r * fly3i;
  2095.             fly1i = fly0i + U1r * fly3r;
  2096.             fly1i = fly1i + U1i * fly3i;
  2097.             fly0r = scale * fly0r - fly1r;
  2098.             fly0i = scale * fly0i - fly1i;
  2099.  
  2100.             fly6r = t0r - U2r * fly4r;
  2101.             fly6r = fly6r - U2i * fly4i;
  2102.             fly6i = t0i + U2i * fly4r;
  2103.             fly6i = fly6i - U2r * fly4i;
  2104.             t0r = scale * t0r - fly6r;
  2105.             t0i = scale * t0i - fly6i;
  2106.  
  2107.             fly3r = fly5r - U2i * fly2r;
  2108.             fly3r = fly3r + U2r * fly2i;
  2109.             fly3i = fly5i - U2r * fly2r;
  2110.             fly3i = fly3i - U2i * fly2i;
  2111.             fly2r = scale * fly5r - fly3r;
  2112.             fly2i = scale * fly5i - fly3i;
  2113.  
  2114.             fly4r = t1r - U3r * fly0r;
  2115.             fly4r = fly4r - U3i * fly0i;
  2116.             fly4i = t1i + U3i * fly0r;
  2117.             fly4i = fly4i - U3r * fly0i;
  2118.             t1r = scale * t1r - fly4r;
  2119.             t1i = scale * t1i - fly4i;
  2120.  
  2121.             fly5r = fly7r + U3i * fly1r;
  2122.             fly5r = fly5r - U3r * fly1i;
  2123.             fly5i = fly7i + U3r * fly1r;
  2124.             fly5i = fly5i + U3i * fly1i;
  2125.             fly7r = scale * fly7r - fly5r;
  2126.             fly7i = scale * fly7i - fly5i;
  2127.  
  2128.             *(fly0P + FlyOffsetA) = fly6r;
  2129.             *(fly0P + FlyOffsetAIm) = fly6i;
  2130.             *(fly0P) = t0r;
  2131.             *(fly0P + 1) = t0i;
  2132.             *(fly2P) = fly3r;
  2133.             *(fly2P + 1) = fly3i;
  2134.             *(fly2P + FlyOffsetA) = fly2r;
  2135.             *(fly2P + FlyOffsetAIm) = fly2i;
  2136.  
  2137.             fly0r = *(fly0P + FlyOffsetB);
  2138.             fly0i = *(fly0P + FlyOffsetBIm);
  2139.  
  2140.             *(fly1P + FlyOffsetA) = fly4r;
  2141.             *(fly1P + FlyOffsetAIm) = fly4i;
  2142.             *(fly1P) = t1r;
  2143.             *(fly1P + 1) = t1i;
  2144.  
  2145.             fly1r = *(fly1P + FlyOffsetB);
  2146.             fly1i = *(fly1P + FlyOffsetBIm);
  2147.  
  2148.             *(fly3P + FlyOffsetA) = fly5r;
  2149.             *(fly3P + FlyOffsetAIm) = fly5i;
  2150.             *(fly3P) = fly7r;
  2151.             *(fly3P + 1) = fly7i;
  2152.  
  2153.             fly0P = (fly0P + FlyOffsetB);
  2154.             fly1P = (fly1P + FlyOffsetB);
  2155.             fly2P = (fly2P + FlyOffsetB);
  2156.             fly3P = (fly3P + FlyOffsetB);
  2157.         };
  2158.  
  2159.         fly2r = *(fly2P);
  2160.         fly2i = *(fly2P + 1);
  2161.         fly3r = *(fly3P);
  2162.         fly3i = *(fly3P + 1);
  2163.         fly4r = *(fly0P + FlyOffsetA);
  2164.         fly4i = *(fly0P + FlyOffsetAIm);
  2165.         fly5r = *(fly1P + FlyOffsetA);
  2166.         fly5i = *(fly1P + FlyOffsetAIm);
  2167.         fly6r = *(fly2P + FlyOffsetA);
  2168.         fly6i = *(fly2P + FlyOffsetAIm);
  2169.         fly7r = *(fly3P + FlyOffsetA);
  2170.         fly7i = *(fly3P + FlyOffsetAIm);
  2171.  
  2172.         t1r = fly0r - U0r * fly1r;
  2173.         t1r = t1r - U0i * fly1i;
  2174.         t1i = fly0i + U0i * fly1r;
  2175.         t1i = t1i - U0r * fly1i;
  2176.         t0r = scale * fly0r - t1r;
  2177.         t0i = scale * fly0i - t1i;
  2178.  
  2179.         fly1r = fly2r - U0r * fly3r;
  2180.         fly1r = fly1r - U0i * fly3i;
  2181.         fly1i = fly2i + U0i * fly3r;
  2182.         fly1i = fly1i - U0r * fly3i;
  2183.         fly2r = scale * fly2r - fly1r;
  2184.         fly2i = scale * fly2i - fly1i;
  2185.  
  2186.         fly0r = fly4r - U0r * fly5r;
  2187.         fly0r = fly0r - U0i * fly5i;
  2188.         fly0i = fly4i + U0i * fly5r;
  2189.         fly0i = fly0i - U0r * fly5i;
  2190.         fly4r = scale * fly4r - fly0r;
  2191.         fly4i = scale * fly4i - fly0i;
  2192.  
  2193.         fly3r = fly6r - U0r * fly7r;
  2194.         fly3r = fly3r - U0i * fly7i;
  2195.         fly3i = fly6i + U0i * fly7r;
  2196.         fly3i = fly3i - U0r * fly7i;
  2197.         fly6r = scale * fly6r - fly3r;
  2198.         fly6i = scale * fly6i - fly3i;
  2199.  
  2200.  
  2201.         fly5r = t0r - U1r * fly2r;
  2202.         fly5r = fly5r - U1i * fly2i;
  2203.         fly5i = t0i + U1i * fly2r;
  2204.         fly5i = fly5i - U1r * fly2i;
  2205.         t0r = scale * t0r - fly5r;
  2206.         t0i = scale * t0i - fly5i;
  2207.  
  2208.         fly7r = t1r + U1i * fly1r;
  2209.         fly7r = fly7r - U1r * fly1i;
  2210.         fly7i = t1i + U1r * fly1r;
  2211.         fly7i = fly7i + U1i * fly1i;
  2212.         t1r = scale * t1r - fly7r;
  2213.         t1i = scale * t1i - fly7i;
  2214.  
  2215.         fly2r = fly4r - U1r * fly6r;
  2216.         fly2r = fly2r - U1i * fly6i;
  2217.         fly2i = fly4i + U1i * fly6r;
  2218.         fly2i = fly2i - U1r * fly6i;
  2219.         fly4r = scale * fly4r - fly2r;
  2220.         fly4i = scale * fly4i - fly2i;
  2221.  
  2222.         fly1r = fly0r + U1i * fly3r;
  2223.         fly1r = fly1r - U1r * fly3i;
  2224.         fly1i = fly0i + U1r * fly3r;
  2225.         fly1i = fly1i + U1i * fly3i;
  2226.         fly0r = scale * fly0r - fly1r;
  2227.         fly0i = scale * fly0i - fly1i;
  2228.  
  2229.         U0i = *(U0iP = (U0iP - NsameU4));
  2230.         U0r = *(U0rP = (U0rP + NsameU4));        
  2231.         U1r = *(U1rP = (U1rP + NsameU2));
  2232.         U1i = *(U1iP = (U1iP - NsameU2));
  2233.         if(stage&1)
  2234.             U0r = -U0r;
  2235.  
  2236.         fly6r = t0r - U2r * fly4r;
  2237.         fly6r = fly6r - U2i * fly4i;
  2238.         fly6i = t0i + U2i * fly4r;
  2239.         fly6i = fly6i - U2r * fly4i;
  2240.         t0r = scale * t0r - fly6r;
  2241.         t0i = scale * t0i - fly6i;
  2242.  
  2243.         fly3r = fly5r - U2i * fly2r;
  2244.         fly3r = fly3r + U2r * fly2i;
  2245.         fly3i = fly5i - U2r * fly2r;
  2246.         fly3i = fly3i - U2i * fly2i;
  2247.         fly2r = scale * fly5r - fly3r;
  2248.         fly2i = scale * fly5i - fly3i;
  2249.  
  2250.         fly4r = t1r - U3r * fly0r;
  2251.         fly4r = fly4r - U3i * fly0i;
  2252.         fly4i = t1i + U3i * fly0r;
  2253.         fly4i = fly4i - U3r * fly0i;
  2254.         t1r = scale * t1r - fly4r;
  2255.         t1i = scale * t1i - fly4i;
  2256.  
  2257.         fly5r = fly7r + U3i * fly1r;
  2258.         fly5r = fly5r - U3r * fly1i;
  2259.         fly5i = fly7i + U3r * fly1r;
  2260.         fly5i = fly5i + U3i * fly1i;
  2261.         fly7r = scale * fly7r - fly5r;
  2262.         fly7i = scale * fly7i - fly5i;
  2263.  
  2264.         *(fly0P + FlyOffsetA) = fly6r;
  2265.         *(fly0P + FlyOffsetAIm) = fly6i;
  2266.         *(fly0P) = t0r;
  2267.         *(fly0P + 1) = t0i;
  2268.  
  2269.         U2r = *(U2rP = (U2rP + NsameU1));
  2270.         U2i = *(U2iP = (U2iP - NsameU1));
  2271.  
  2272.         *(fly2P) = fly3r;
  2273.         *(fly2P + 1) = fly3i;
  2274.         *(fly2P + FlyOffsetA) = fly2r;
  2275.         *(fly2P + FlyOffsetAIm) = fly2i;
  2276.         *(fly1P + FlyOffsetA) = fly4r;
  2277.         *(fly1P + FlyOffsetAIm) = fly4i;
  2278.         *(fly1P) = t1r;
  2279.         *(fly1P + 1) = t1i;
  2280.  
  2281.         U3r = *( U2rP + U3offset);
  2282.         U3i = *( U2iP - U3offset);
  2283.  
  2284.         *(fly3P + FlyOffsetA) = fly5r;
  2285.         *(fly3P + FlyOffsetAIm) = fly5i;
  2286.         *(fly3P) = fly7r;
  2287.         *(fly3P + 1) = fly7i;
  2288.  
  2289.         IOP = IOP + 2;
  2290.         fly0P = IOP;
  2291.         fly1P = (IOP+Flyinc);
  2292.         fly2P = (fly1P+Flyinc);
  2293.         fly3P = (fly2P+Flyinc);
  2294.     };
  2295.     NsameU4 = -NsameU4;
  2296.  
  2297.     if(stage&1){
  2298.         LoopN >>= 3;
  2299.         NsameU1 >>= 3;
  2300.         NsameU2 >>= 3;
  2301.         NsameU4 >>= 3;
  2302.         NdiffU <<= 3;
  2303.         Flyinc = Flyinc << 3;
  2304.         FlyOffsetA <<= 3;
  2305.         FlyOffsetB <<= 3;
  2306.         FlyOffsetAIm =  FlyOffsetA + 1;
  2307.         FlyOffsetBIm =  FlyOffsetB + 1;
  2308.     }
  2309. }
  2310.  
  2311. /*    Finish RFFT        */
  2312. M=M+1;
  2313.  
  2314. FlyOffsetA = Ntbl[M] * 2/4;
  2315. FlyOffsetAIm = Ntbl[M] * 2/4 + 1;
  2316.  
  2317. IOP = ioptr;
  2318.  
  2319. fly0P = (IOP + Ntbl[M]*2/4);
  2320. fly1P = (IOP + Ntbl[M]*2/8);
  2321.  
  2322. U0rP = &Utbl[Ntbl[M-3]];
  2323.  
  2324. U0r =  *U0rP,
  2325.  
  2326. fly0r = *(IOP);
  2327. fly0i = *(IOP + 1);
  2328. fly1r = *(fly0P);
  2329. fly1i = *(fly0P + 1);
  2330. fly2r = *(fly1P);
  2331. fly2i = *(fly1P + 1);
  2332. fly3r = *(fly1P + FlyOffsetA);
  2333. fly3i = *(fly1P + FlyOffsetAIm);
  2334.  
  2335.     t0r = scale * fly0r + scale * fly0i;
  2336.     t0i = 0.0;
  2337.     t1r = scale * fly1r;
  2338.     t1i = -scale * fly1i;
  2339.  
  2340.     fly0r = fly2r + fly3r;
  2341.     fly0i = fly2i - fly3i;
  2342.     fly1r = fly2i + fly3i;
  2343.     fly1i = fly3r - fly2r;
  2344.  
  2345.     fly2r = fly0r + U0r * fly1r;
  2346.     fly2r = fly2r + U0r * fly1i;
  2347.     fly2i = fly0i - U0r * fly1r;
  2348.     fly2i = fly2i + U0r * fly1i;
  2349.     fly3r = scale * fly0r - fly2r;
  2350.     fly3i = fly2i - scale * fly0i;
  2351.  
  2352. *(IOP) = t0r;
  2353. *(IOP + 1) = t0i;
  2354. *(fly0P) = t1r;
  2355. *(fly0P + 1) = t1i;
  2356. *(fly1P) = fly2r;
  2357. *(fly1P + 1) = fly2i;
  2358. *(fly1P + FlyOffsetA) = fly3r;
  2359. *(fly1P + FlyOffsetAIm) = fly3i;
  2360.  
  2361. U0rP = &Utbl[1];
  2362. U0iP = &Utbl[Ntbl[M-2]-1];
  2363.  
  2364. U0r =  *U0rP,
  2365. U0i =  *U0iP;
  2366.     
  2367. fly0P = (IOP + 2);
  2368. fly1P = (IOP + (Ntbl[M-2]-1)*2);
  2369.  
  2370.                 /* Butterflys */
  2371.                 /*
  2372.         f0    -    t0    -    -    f2
  2373.         f1    -    t1    -U0    -    f3
  2374.         f2    -    f0    -    -    t0
  2375.         f3    -    f1    -U0a-    t1
  2376.                 */
  2377.  
  2378. for (diffUcnt = Ntbl[M-3]-1; diffUcnt > 0 ; diffUcnt--){
  2379.  
  2380.     fly0r = *(fly0P);
  2381.     fly0i = *(fly0P + 1);
  2382.     fly1r = *(fly1P + FlyOffsetA);
  2383.     fly1i = *(fly1P + FlyOffsetAIm);
  2384.     fly2r = *(fly1P);
  2385.     fly2i = *(fly1P + 1);
  2386.     fly3r = *(fly0P + FlyOffsetA);
  2387.     fly3i = *(fly0P + FlyOffsetAIm);
  2388.  
  2389.     t0r = fly0r + fly1r;
  2390.     t0i = fly0i - fly1i;
  2391.     t1r = fly0i + fly1i;
  2392.     t1i = fly1r - fly0r;
  2393.  
  2394.     fly0r = fly2r + fly3r;
  2395.     fly0i = fly2i - fly3i;
  2396.     fly1r = fly2i + fly3i;
  2397.     fly1i = fly3r - fly2r;
  2398.  
  2399.     fly2r = t0r + U0r * t1r;
  2400.     fly2r = fly2r + U0i * t1i;
  2401.     fly2i = t0i - U0i * t1r;
  2402.     fly2i = fly2i + U0r * t1i;
  2403.     fly3r = scale * t0r - fly2r;
  2404.     fly3i = fly2i - scale * t0i;
  2405.  
  2406.     t0r = fly0r + U0i * fly1r;
  2407.     t0r = t0r + U0r * fly1i;
  2408.     t0i = fly0i - U0r * fly1r;
  2409.     t0i = t0i + U0i * fly1i;
  2410.     t1r = scale * fly0r - t0r;
  2411.     t1i = t0i - scale * fly0i;
  2412.  
  2413.     *(fly0P) = fly2r;
  2414.     *(fly0P + 1) = fly2i;
  2415.     *(fly1P + FlyOffsetA) = fly3r;
  2416.     *(fly1P + FlyOffsetAIm) = fly3i;
  2417.  
  2418.     U0r = *(U0rP = (U0rP + 1));        
  2419.     U0i = *(U0iP = (U0iP - 1));
  2420.  
  2421.     *(fly1P) = t0r;
  2422.     *(fly1P + 1) = t0i;
  2423.     *(fly0P + FlyOffsetA) = t1r;
  2424.     *(fly0P + FlyOffsetAIm) = t1i;
  2425.  
  2426.     fly0P += 2;
  2427.     fly1P -= 2;
  2428. };
  2429.  
  2430. ioptr += Ntbl[M];
  2431. }
  2432. }
  2433.